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, ABSTRACT 

■ We present predictions for the reionization of the intergalactic medium (IGM) by 
I stars in high-redshift galaxies, based on a semi-analytic model of galaxy formation. 

■ We calculate ionizing luminosities of galaxies, including the effects of absorption by 
^ ' interstellar gas and dust on the escape fraction /esc, and follow the propagation of 

I the ionization fronts around each galaxy in order to calculate the filling factor of 

■ ionized hydrogen in the IGM. For a ACDM cosmology, with parameters of the galaxy 
' formation model chosen to match observations of present-day galaxies, and a physical 

calculation of the escape fraction, we find that the hydrogen in the IGM will be 
reionized at redshift z = 6.1 if the IGM has uniform density, but only by z = 4.5 if the 
IGM is clumped. If instead we assume a constant escape fraction of 20% for all galaxies, 
then we find reionization at z = 9.0 and z — 7.8 for the same two assumptions about 
i-S^ ' IGM clumping. We combine our semi-analytic model with an N-body simulation of the 

distribution of dark matter in the universe in order to calculate the evolution of the 
Q . spatial and velocity distribution of the ionized gas in the IGM, and use this to calculate 

$-H ' the secondary temperature anisotropics induced in the cosmic microwave background 

I (CMB) by scattering off free electrons. The models predict a spectrum of secondary 

■ anisotropics covering a broad range of angular scales, with fractional temperature 
^ ' fluctuations ^ 10~^ — 10~^ on arcminute scales. The amplitude depends strongly on 
i-H . the total baryon density, and less sensitively on the escape fraction /osc- The amplitude 
^ ' also depends somewhat on the geometry of reionization, with models in which the 
H I regions of highest gas density are reionized first giving larger CMB fluctuations than 
. . 1 the case where galaxies ionize surrounding spherical regions, and models where low 

density regions reionize first giving the smallest fluctuations. Measurement of these 
anisotropics can therefore put important constraints on the reionization process, in 
particular, the redshift evolution of the filling factor, and should be a primary objective 
of a next generation submillimeter telescope such as the Atacama Large Millimeter 
Array. 

Key words: cosmology: theory- dark matter- large scale structure of Universe- in- 
tergalactic medium 



1 INTRODUCTION 

The Gunn-Peterson (GP) effect (Gunn & Peterson 1965) 
strongly indicates that the smoothly distributed hydrogen 
in the intergalactic medium (IGM) is already highly ion- 
ized by z = 5 (Schneider, Schmidt & Gunn 1991; Lanzetta, 
Wolfe & Turnshek 1996). Barring the possibility of colli- 
sional reionization (e.g. Giroux & Shapiro 1994), the GP 
effect implies the presence of very luminous ionizing sources 



at high redshifts capable of producing enough Lyman con- 
tinuum (Lyc) photons to cause photoionization of hydrogen 
by z ^5. The two possible sources of these ionizing photons 
are QSOs and high mass stars. 

Models in which QSOs dominate the production of 
ionizing photons may be able to meet the GP constraint 
(Miralda-Escude & Ostriker 1990). However, such models 
are strongly constrained by the observed drop in the abun- 
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dance of bright QSOs above 2 « 3 (Hartwick & Schade 1990; 
Warren, Hewett & Osmer 1994; Kennefick, Djorgovski & de 
Carvalho 1995; Schmidt, Schneider & Gunn 1995). Further- 
more, Madau, Haardt & Rees (1999) note that a model in 
which faint QSOs provide all the required ionizing luminos- 
ity can be ruled out on the basis of the number of faint QSOs 
seen in the HDF. 

There is, however, growing evidence for the presence of 
bright galaxies at redshifts as high as z ~ 5 (Spinrad 1998), 
and perhaps even higher (Yahil, Lanzetta & Ferna'ndez-Soto 
1998). Thus the other natural candidate sources of ioniz- 
ing photons are young, high mass stars forming in galaxies 
at redshifts greater than 5 (e.g. Couchman & Rees 1986; 
Haiman&Loeb 1996; Ciardi et al. 2000) . Madau, Haardt & 
Rees (1999) note that at 2; « 3 stars in Lyman-break galax- 
ies will emit more ionizing photons into the IGM than QSOs 
if more than 30% of such photons can escape from their host 
galaxy. Whilst such high escape fractions may not be real- 
istic (e.g. local starbursts show escape fractions of only a 
few percent, Leitherer et al. 1995), this does demonstrate 
that high-redshift galaxies could provide a significant con- 
tribution to (or perhaps even dominate) the production of 
ionizing photons. In this work we will restrict our attention 
to ionizing photons produced by stars, deferring considera- 
tion of the QSO contribution to a later paper. 

According to the hierarchical structure formation sce- 
nario (e.g. Peebles 1980) perturbations in the gravitation- 
ally dominant and dissipationless dark matter grow, by grav- 
itational instability, into virialised clumps, or halos. Galax- 
ies, and later stars, then form by the cooling and conden- 
sation of gas inside these halos (e.g. White & Rees 1978; 
White & Frenk 1991). Dark matter halos continually grow 
by merging with other halos (e.g. Bower 1991; Bond, Cole, 
Efstathiou & Kaiser 1991). In the context of this hierar- 
chical scenario, we present a realistic scheme for studying 
the reionization of the universe by ionizing photons emitted 
from massive stars. We focus on the photoionization of the 
hydrogen component of the IGM. To predict the time depen- 
dent luminosity in Lyc photons we use a semi-analytic model 
of galaxy formation (e.g. Kauffmann et al. 1993; Cole et al. 
1994; Somerville & Primack 1999). In particular, we use the 
semi-analytic model of Cole et al. (2000), modified to take 
into account Compton cooling by cosmic microwave back- 
ground (CMB) photons, to model the properties of galaxies 
living in dark matter halos spanning a wide range of masses. 

We then estimate the fraction of the ionizing photons 
which manage to escape each galaxy, and therefore con- 
tribute to the photoionization of the intergalactic Hi. The 
fraction of ionizing photons escaping is determined on a 
galaxy- by-galaxy basis, using physically motivated models. 
Assuming spherical symmetry, we follow the propagation of 
the ionization front around each halo to compute the fill- 
ing factor of intergalactic Hll regions, including the efi'ects 
of clumping in the IGM. Finally, using several alternative 
models for the spatial distribution of ionized regions within a 
high resolution N-body simulation of the dark matter distri- 
bution, we estimate the anisotropies imprinted on the CMB 
by the patchy reionization process, due to the correlations 
in the ionized gas distribution and velocities (Sunyaev & 
Zel'dovich 1980; Vishniac 1987). In previous models many 
simplifications were made in computing both the spatial dis- 
tribution of ionized regions and the two-point correlations 



of gas density and velocity in those regions (Aghanim et al. 
1996; Jaffe & Kamionkowski 1998; Gruzinov & Hu 1998; 
Knox, Scoccimarro & Dodelson 1998; Peebles & Juskiewicz 
1998; Haiman & Knox 1999). Our calculations represent a 
significant improvement over these models as we are able to 
calculate the two-point correlations between gas density and 
velocity in ionized regions directly from an N-body simula- 
tion. 

The rest of this paper is arranged as follows. In §^ we 
outline the features of the semi-analytic model relevant to 
galaxy formation at high redshifts. In §^ we describe how 
we calculate the fraction of ionizing photons escaping from 
galaxies, and observational constraints on the ionizing lumi- 
nosities and escape fraction at low and high redshift from 
Ha luminosities and Hi masses and column densities. In §^ 
we describe how we calculate the filling factor of photoion- 
ized gas in the IGM, including the effects of clumping of this 
gas. We then present our predictions for reionization, includ- 
ing the effects on the reionization redshift of using different 
assumptions about escape fractions and clumping factors. In 
^ we examine the robustness of our results to changes in 
the other parameters of the semi-analytic galaxy formation 
model. In ^ we describe how the semi-analytic models are 
combined with N-body simulations to calculate the spatial 
distribution of the photoionized IGM. We then calculate the 
spectrum of anisotropies introduced into the CMB by this 
ionized gas. Finally, in we summarize our results and ex- 
amine their consequences. 



2 THE SEMI-ANALYTIC MODEL OF 
GALAXY FORMATION 

To determine the luminosity in ionizing Lyc photons pro- 
duced by the galaxy population, we use the semi-analytic 
model of galaxy formation developed by Cole et al. (2000). 
This model predicts the properties of galaxies residing 
within dark matter halos of different masses. This is achieved 
by relating, in a self-consistent way, the physical processes 
of gas cooling, star formation, and supernovae feedback to a 
halo's merger history, which is calculated using the extended 
Press-Schechter theory. The parameters of this model are 
constrained by a set of observations of galaxies in the local 
Universe, including the B and K-band luminosity functions, 
the I-band TuUy-Fisher relation, the mixture of morphologi- 
cal types and the distribution of disk scale lengths (see Cole 
et al. (2000) and references therein for a thorough discus- 
sion of the observational constraints). Once the model has 
been constrained in this way it is able to make predictions 
concerning the clustering of galaxies (Benson et al. 2000) 
and the properties of galaxies at higher redshifts. For ref- 
erence, the parameters of our standard model are given in 
Table |^. Definitions of the semi-analytic model parameters 
can be found in Cole et al. (2000). 

As we are employing the semi-analytic model at much 
higher redshifts than we have previously attempted, we will 
investigate the effects on our results of changing key model 
parameters. Of particular interest will be the prescription 
for feedback from supernovae and stellar winds. The model 
assumes that a mass /3AM of gas is reheated by supernovae 
and ejected from the disk for each mass AM of stars formed. 
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Table 1. The parameters of our standard model. 



Parameter 


Value 


Cosmology 


Qo 


0.3 


Ao 


0.7 


Hq 


70 km/s/Mpc 


OS 


0.90 


r 


0.21 




0.02 



Gas cooling 



Gas profile^ 

^core 

Recooling* 



CDC 

0.33rNFw 
not allowed 



Star formation and feedback 





-1.5 




0.01 


ahot 


2.0 




150 km s-l 


Feedback** standard 


Stellar populations 


IMF 


Kennicutt (1983) 


P 


0.02 


R 


0.31 


T 


1.53 


Mergers and bursts 


fdf 


1.0 


/ellip 


0.3 


Starbursts included 


/dyn 


1.0 


Ionizing luminosities 




10.0 X 10^° photons/s 




0.1 



^ Gas profiles that we consider are CDC (for which the gas den- 
sity is p(r) oc (r^ +'"?orc)~^ — i-^- isothermal profile but with 
a constant density core), and SIS (an isothermal profile with no 
core). 

* Gas ejected from galaxies by supernovae is allowed to cool again 
in the same halo if 'Recooling' is allowed. Otherwise this gas can 
only cool again once it enters a newly formed halo. 

** 'Standard' feedback is the form specified in eqn. (|l|). The alter- 
native is 'modified' feedback, which is the form specified in eqn. 



t The upper cut off the in luminosity function of O B as sociations 
assumed in the DS94 model (see Appendix ^ and j ^.2| ). 

t This is the ratio of disk vertical and radial scale lengths used 
in the DS94 and DSG N98 models for the gas escape fraction (see 
Appendix ^ and [3.2). 



The quantity /3 is allowed to be a function of the galaxy 
properties, and is parameterised as 



/3= (VdiBk/Vhot)" 



(1) 



where Vdisk is the circular velocity of the galaxy disk and 
Vhot and Qhot are adjustable parameters of the model. Cole 
et al. (2000) show that Qhot and Viiot are well constrained 
by the shape of the B-band luminosity function and the 
TuUy-Fisher relation at 2 = 0. However, since there is very 
little time available for star formation at the high redshifts 
which we consider, it is possible that these parameters, or 
indeed the form of the parameterisation in eqn. (|l|), could 
be changed at high redshift without significantly affecting 
the model predictions at z = 0. In ^ we will therefore ex- 
periment with different values of these parameters and will 
also consider a modified functional form for /?. 

Two other key model inputs are the baryon density pa- 
rameter, f2b, and the stellar initial mass function (IMF). The 
value of f2b determines cooling rates (and so star formation 
rates) in our model halos. The shape of the IMF determines 
the number of high mass stars which produce the ionizing 
photons. For fib our standard value is 0.02, which is consis- 
tent with the estimate Walker et al. (1991) and which allows 
a good match to the bright end of the observed B-band lumi- 
nosity function. We will also consider an alternative value of 
Qh = 0.04, which is in better agreement with estimates from 
the D/H ratio in QSO absorption line systems (Schramm & 
Turner 1998; Buries & Tytler 1998). For the IMF, we adopt 
as our standard choice the IMF of Kennicutt (1983), which 
is close to the "best" IMF proposed by Scalo (1998) on 
the basis of observations in the Solar neighbourhood and in 
nearby galaxies. We consider the effects of changing both 57b 
and the IMF in sH. 



2.1 Gas Cooling 

The standard semi-analytic model of Cole et al. (2000) al- 
lows hot halo gas to cool only via collisional radiative pro- 
cesses. At high redshift, Compton cooling due to free elec- 
trons in the hot plasma scattering off CMB photons becomes 
important. The Compton cooling timescale is given by (Pee- 
bles 1968) 



^Co 



1161.3(1 + x~ 



(l+^)4(l_rCMB(i_^^)/y^) 



Gyr, 



(2) 



where a;hot is the ionized fraction of the hot halo gas, T^^^ 
is the temperature of the CMB at the present day and To is 
the temperature of electrons in the hot halo gas, which we set 
equal to the virial temperature of the halo. At high redshifts 
this cooling time becomes shorter than the Hubble time and 
so Compton cooling may be effective at these redshifts. To 
implement eqn. (^) in the semi-analytic model, we assume 
that the shock-heated halo gas is in collisional ionization 
equilibrium, and use values for a::hot which we interpolate 
from the tabulated values given by Sutherland & Dopita 
(1993). In halos with virial temperatures less than around 
lO'^K collisional ionization is ineffective, and so the ionized 
fraction in the halo gas will equal the residual ionization 
fraction left over from recombination. However, since this 
fraction is small, we will simply assume in this paper that 
cooling in halos below lO^K is negligible. 
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It should be noted that, unhke the radiative cooling 
time, the Compton cooling time is independent of the gas 
density and depends only very weakly on the gas temper- 
ature. Whereas with coUisional radiative cooling a cooling 
radius, within which the cooling time is less than the age 
of the halo, propagates through the halo as more and more 
gas cools, with Compton cooling the entire halo cools at the 
same rate. The amount of gas which can reach the centre of 
the halo is then controlled by the free-fall timescale in the 
halo. 

Including Compton cooling in our model turns out to 
make little difference to the results. For example, the total 
mass of stars formed in the universe as a function of redshift 
differs by less than 5% for z < 20 between models with and 
without Compton cooling. At higher redshifts the differences 
can become as large as 30% for small intervals of redshift. 
For example, if a massive halo cools via Compton cooling it 
will rapidly produce many stars — without Compton cooling 
it will still form these stars, but not until slightly later when 
coUisional radiative cooling takes effect. However, the mass 
of stars formed at these high redshifts is tiny, and so any 
differences become entirely negligible at lower redshifts when 
many more stars have formed. Although at high redshift 
the Compton cooling time is shorter than the age of the 
Universe, halos are merging at a high rate and so their gas is 
being repeatedly shock-heated by successive mergers which, 
we assume, heat the gas to the virial temperature of the 
halo. We find that for the majority of halos the time between 
successive major mergers (defined as the time for a doubling 
in mass of the halo) is less than the Compton cooling time 
at the redshifts considered here. Therefore, Compton cooling 
will be ineffective in these halos. In the few cases where the 
halo does survive long enough that Compton cooling could 
be important, we find that the coUisional radiative cooling 
time at the virial radius of the halo is often shorter than 
the Compton cooling time, in which case all of the gas in 
the halo will cool whether or not we include the effects of 
Compton cooling. Nevertheless, Compton cooling is included 
in all models considered. 

In this work we ignore cooling due to molecular hydro- 
gen (H2). Although molecular hydrogen allows cooling to 
occur in gas below lO'^K, it is easily dissociated by photons 
from stars that form from the cooling gas. Previous studies 
that have included cooling due to H2 typically find that it 
is completely dissociated at very high redshifts. For exam- 
ple, Ciardi et al. (2000) find that molecular hydrogen is 
fully dissociated by z « 25. Objects formed by II2 cooling 
are therefore not expected to contribute significantly to the 
reionization of the IGM. 



2.2 Fraction of Gas in the IGM 

At any given redshift, some fraction of the gas in the Uni- 
verse will have become coUisionally ionized in dark matter 
halos and some fraction will have cooled to become part of 
a galaxy. Within the context of our semi-analytic model, we 
define the IGM as all gas which has not been coUisionally 
ionized inside dark matter halos and which has not become 
part of a galaxy (note that we are here only interested in ion- 
ization of hydrogen) . It is this gas which must be photoion- 
ized if the Gunn-Peterson constraint is to be satisfied. The 
fraction of the total baryon content of the universe which is 




0.8 - 



0.6 



0.4 



0.2 



Figure 1. The fraction of baryons remaining in the IGM as a 
function of redshift. Other baryons have either been coUisionally 
ionized in dark matter halos or have cooled to become part of a 
galaxy. 



in the IGM, /igm, can be estimated by integrating over the 
mass function of dark matter halos, as follows 



/igm(z) = 1- 



dn dMhaio 



dMh 



/galaxy (2), (3) 



where A/gas(Mhaio, 2) is the mean mass of diffuse gas in ha- 
los of mass Mhaio, xniMhaio, z) is the fraction of hydrogen 
which is coUisionally ionized at the halo virial temperature 
(which we take from the calculations of Sutherland & Dopita 
1993), dn/dAfhaio(Mhaio, z) is the comoving number density 
of halos (which we approximate by the Press-Schechter mass 
function), pc is the critical density of the Universe at z = 0, 
and /galaxy IS the fraction of the total baryonic mass in the 
Universe which has been incorporated into galaxies. The 
quantities Mgas and /galaxy can be readily calculated from 
our model of galaxy formation. Fig. ^ shows the evolution 
of ,fiGM with redshift. 



2.3 Observational Constraints 

The semi-analytic model provides the spectral energy dis- 
tribution (SED) of each galaxy, from which we can deter- 
mine the ionizing luminosity of that galaxy. Summing the 
contributions from all galaxies in a given halo yields the to- 
tal ionizing luminosity produced in that halo. Cojazzi et al. 
(2000), using the model of Haiman & Loeb (1996), demon- 
strated that a higher reionization redshift could be obtained 
if zero-metallicity stars were responsible for reionization, as 
these produce a greater number of ionizing photons than 
low (i.e. 10~*) metallicity stars. In our model the very first 
stars have zero metallicity, but as we include chemical evo- 
lution only a very small fraction of stars have metallicities 
below 10"*. This is consistent with the results of Tumlinson 
& ShuU (2000) who argue that the epoch of metal-free star 
formation must end before 2; = 3, as the enhanced emis- 
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Figure 2. The comoving number density of galaxies brighter than 
a given ionizing luminosity function of redshift. The 

number is plotted for nion = 10^ (solid line), 10^ (dotted line) 
and 10* (dashed line) in units of 10^" photons/s. 

sion shortwards of 228A from such stars is inconsistent with 
observations of Hell opacity in the IGM at that redshift. 
Therefore we cannot appeal to such zero-metaUicity stars to 
increase the redshift of reionization in our modeL 

As a result of absorption by neutral hydrogen close to 
the emitting stars and extinction caused by dust, only a 
small fraction of the ionizing radiation emitted by the stars 
escapes from each galaxy (Leitherer et al. 1995; Hurwitz, 
Jehnsky & Dixon 1997; Kunth et al. 1998). We therefore 
estimate, within the context of the semi-analytic model, the 
fraction, /esc, of ionizing photons which escape the galaxy 
to become available for the photoionization of the Hi in the 
IGM. The calculation of /esc is discussed in 

In Fig. P we show the redshift evolution of the comov- 
ing number density of galaxies with ionizing luminosity hion 
larger than 10^, lO'^ and 10* in units of 10^" photons per 
second. These are the unattenuated luminosities produced 
by massive stars in the galaxies. The abundances of sources 
of given luminosity rises sharply up to 2: = 2 — 4 (the exact 
position of the peak depending on luminosity) as more and 
more dark matter halos form that are capable of hosting 
bright galaxies. After 2 = 2 — 4 abundances quickly drop to- 
wards 2 = as the amount of gas available for star formation 
declines. 

The escape fractions in our model will be determined 
by the mass and radial scale length of the Hi gas in 
galactic disks. It is therefore important to test that our 
model produces galaxies with reasonable distributions of 
Hi mass and disk scale length. Cole et al. (2000) have 
shown that our semi-analytic model produces distributions 
of I-band disk scale lengths in good agreement with the 
2 = data of de Jong & Lacey (1999). In Fig. || we 
compare our model with observations of damped Lyman- 
Q systems (DLAS) over a range of redshifts and with the 
Hi mass function at 2 = 0.0. Under the assumption that 




log(N„,/cm-2) log(M„,/h-2M,J 



Figure 3. Panels (a), (b) and (c) show the distribution of DLAS 
at 2 = 1, 2 and 3 respectively as a function of their Hi column 
density. Solid lines indicate the distribution determined from our 
model. Points with error bars are from Lanzetta, Wolfe & Turn- 
shek (1996). Panel (d) shows the Hi mass function of galaxies at 
z = 0.0. The solid line is the mass function determined from our 
model. Points with errorbars are from Zwaan ct al. (1997). 



DLAS are caused by neutral gas in galactic disks, we com- 
pute the DLAS column density distribution in our model, 
/dlas, defined such that /dlas(A'hi, i)dA'HidA is the mean 
number of DLAS at cosmic time t with column densities 
in the range A'^hi to A''hi -I- dA^Hi and absorption distance 
X{z) = f [(1 + 2)^/^ - 1] in the interval dX along a line 
of sight (Lanzetta, Wolfe & Turnshek 1996). Our model 
is in reasonable agreement with the distribution of DLAS 
column densities observed by Lanzetta, Wolfe & Turnshek 
(1996), indicating that both the mass of Hi and its radial 
scale length in our model galaxies are realistic. The 2 = 
Hi mass function from our model is also in reasonable agree- 
ment with the data of Zwaan et al. (1997), although it 
does overpredict the abundance of low Hi mass galaxies. 
Our model predictions assume that all of the hydrogen in 
galactic disks is in the form of Hi. In practice, some of the 
hydrogen in disks will be in the form of molecules (H2) or 
ionized gas (Hll), so this over-estimates the Hi masses and 
column densities. 

A significant contribution to the ionizing luminosity 
comes from very low mass halos. We therefore ensure that 
we resolve all halos which have a virial temperature > lO^K 
up to 2 = 50, i.e. all halos down to a mass of 5 x 10^/i~^Mq. 
Below this temperature cooling becomes inefficient (since we 
are ignoring cooling by molecular hydrogen, and the Comp- 
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ton cooling from the residual free electrons left over after 
recombination) and so galaxy formation ceases. 

The requirement that 5 x 10^ Mq halos be resolved 
sets an upper limit on the mass of halo that we can simulate 
due to computer memory limits, since the lower the mass of 
halo that is resolved, the more progenitors a halo of given 
mass will have. At z = 0, the most massive halos that we are 
able to simulate make a significant contribution to the total 
filling factor and ionizing luminosity. However, for z ~ 2 
the most massive halos simulated contribute only 1% of the 
total number of escaping ionizing photons, and this fraction 
drops extremely rapidly as we look to even higher redshifts. 
Therefore, at the high redshifts {z ^3) we will be interested 
in, ignoring higher mass halos makes no significant difference 
to our results. 

We note that once any halo has begun to ionize the 
surrounding IGM, it could potentially influence the pro- 
cess of galaxy formation in nearby halos. Ionizing photons 
from the first halo will act to heat the gas in nearby ha- 
los, thereby reducing the effective cooling rate (Efstathiou 
1992; Thoul & Weinberg 1996). Since prior to fuU reion- 
ization each halo will see only the fiux of ionizing photons 
from nearby sources, a detailed accounting of this radiative 
feedback requires a treatment of the radiative transfer of 
the ionizing radiation through the IGM. This is beyond the 
scope of the present work. Such radiative feedback is ex- 
pected to be very efficient at dissociating molecular hydro- 
gen, with Ciardi et al. (2000) finding that H2 is completely 
dissociated by z ~ 25. Radiative feedback will also inhibit 
galaxy formation both by reducing the amount of gas that 
accretes into low mass halos (Gnedin 2000) and by reduc- 
ing the cooling rate of gas within halos. Thoul & Weinberg 
(1996) show that radiative feedback may be effective in in- 
hibiting galaxy formation in halos with circular velocities of 
50 km/s or less. In our model, the ionizing luminosity be- 
comes dominated by galaxies in halos with circular velocities 
greater than 50 km/s at redshifts below z « 10. At higher 
redshifts we may therefore be overestimating the total ion- 
izing luminosity produced by galaxies, but this should not 
significantly affect the reionization redshift. 



3 THE ESCAPE FRACTION OF IONIZING 
PHOTONS 

3.1 Global constraints at low redshift 

Gas and dust inside galaxies can readily absorb ionizing pho- 
tons and re-emit the energy at longer wavelengths. Therefore 
the amount and distribution of these components are the 
main factors that determine /osc. The model of galaxy for- 
mation explicitly provides the mass and metallicity of cold 
gas present in each galaxy disk and the half-mass radius of 
that disk, all as functions of time. The mass of dust is as- 
sumed to be proportional to the mass of cold gas and to its 
metallicity. We split the escape fraction into contributions 
from gas, /osc, gas, and dust, /esc, dust, such that the total es- 
caping fraction is given by /esc = /esc, gas /esc, dust. 

Cole et al. (2000) describe in detail how the effects 
of dust are included in their model of galaxy formation. 
This modelling, which uses the calculations of Ferrara et 
al. (1999), is much more realistic than has been previously 



included in semi-analytic galaxy formation models, as it in- 
cludes a fully 3D (though axisymmetric) dust distribution, 
and the dust optical depths are calculated for each galaxy in- 
dividually. In this model, stars are assumed to be distributed 
in a bulge and in an exponential disk with a vertical scale 
height equal to 0.0875 times the radial scale length (this 
ratio was adopted by Ferrara et al. to match the observed 
values for the old disk population of galaxies like the Milky 
Way). The dust is assumed to be distributed in the same 
way as the disk stars. The models give the attenuation of 
the ionizing radiation as a function of the inclination angle 
at which a galaxy is viewed, and we average this over angle 
to find the mean dust extinction for each galaxy. The dust 
attenuations do not include the effects of clumping of the 
stars or dust, and also assume that the ionizing stars have 
the same vertical distribution as the dust. With these two 
caveats in mind, the dust extinctions we apply should only 
be considered as approximate. 

Some of the emitted Lyc photons are absorbed by neu- 
tral hydrogen close to the emitting star, thereby causing Ha 
line emission from the galaxy. Therefore, the Ha luminosity 
function is sensitive to the fraction, /esc, gas, of the ionizing 
photons which manage to escape through the gas. We will 
require our models to reproduce the observed Ha luminosity 
function and luminosity density. 

In Figs. ^ and ^ we compare the Ha emission line prop- 
erties of galaxies in our model with observational data at 
low redshift. The observed values are already corrected for 
dust extinction, so we compare them with the theoretical 
values before dust attenuation. In order to calculate these 
properties accurately, we simulate halos of mass up to and 
including 10^^ Mq. In calculating the Ha line luminos- 
ity of each galaxy, we assume that a fraction 1 — /esc, gas of 
the Lyc photons are absorbed by hydrogen atoms, producing 
Ha photons according to case B recombination. The remain- 
ing Lyc photons escape, after being further attenuated by 
dust. The figures show results for /esc, gas = 0, 0.05 and 0.2, 
which roughly brackets the likely range of values for typical 
disk galaxies at the present day, as we discuss below. Both 
the predicted Ha luminosity function and luminosity density 
are in reasonable agreement with the observations, demon- 
strating that our models produce galaxies with realistic total 
ionizing luminosities (before attenuation by gas and dust). 
In principle, these observational comparisons provide a con- 
straint on the value of /esc, gas, if the other parameters in the 
semi-analytical model are assumed to be known. However, 
in practice it is not possible to reliably distinguish between 
/osc, gas = 0.2 and /esc, gas = Or 0.05, given the uncertainties 
in the observational data. The observational results depend 
on the dust correction factors applied, and there is also some 
uncertainty in the ionizing luminosities predicted by stellar 
population synthesis models for a given IMF. With these 
caveats in mind it would seem that mean escape fractions 
anywhere between zero and 20-30% are acceptable. Obser- 
vations of starburst galaxies in the nearby universe suggest 
that the escape fraction is actually less than 3% for such 
galaxies (Leitherer et al. 1995), but starbursts are known 
to have very high column densities of gas and dust, and 
so the escape fraction in normal galaxies can probably be 
significantly higher (e.g. Kennicutt 1998). 
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Figure 4. The Hq luminosity function at z = 0. Points with er- 
ror bars are observational data from Gallego et al. (1995) (filled 
squares) and Sullivan et al. (2000) (open circles; note that the 
median redshift of this survey is (z) ?s 0.15). The observed Ha 
luminosities are corrected for dust extinction. The solid line is 
the luminosity function from our model assuming that no ioniz- 
ing photons escape (/esc, gas = 0). The dashed line is the same 
function for /osc.gas = 0.05 and the dotted line for /esc, gas = 0.2. 
The Ha luminosities from the model are the values unattcnuated 
by dust. 



3.2 The dependence of fesc.gas on redshift and on 
halo mass 

So far, we have assumed that /esc,gas is a global constant, 
varying neither with galaxy properties nor redshift. The de- 
tails of the physical processes which determine /esc, gas are 
uncertain, but a constant /esc, gas seems unrealistic, as the 
properties of the emitting galaxies depend strongly upon 
both redshift and the mass of the halos in which they live. 
Given the complexity of this problem, here we merely aim 
at establishing the general trend of how /esc, gas may vary 
with halo mass and redshift. 

We will consider three models for /esc, gas. In the first 
model, /esc, gas is assumed to be a universal constant (this 
will be referred to as the "fixed model" ) . In the second and 
third models /esc, gas is evaluated for each galaxy, based on 
its physical properties. These two models are described next. 

Our first physical /esc, gas model is based on the ap- 
proach of Dove & ShuU (1994), hereafter DS94 who de- 
rived an analytic expression for /esc, gas. In their model, Lyc 
photons are emitted by OB associations in a galactic disk 
and escape by ionizing "Hll chimneys" in the Hi layer. The 
fraction of photons escaping a disk of given size and gas 
content can then be calculated. Whilst the original DS94 
model assumes that OB associations all lie in the mid-plane 
of the galaxy disk, we have also considered the case where 
OB associations are distributed vertically like the gas in the 
disk. 

The Hll chimney model of DS94 does not include the 
effects of finite lifetimes of the OB associations, or of dy- 



Figure 5. The Ha luminosity density of the Universe as a func- 
tion of redshift. Points with errorbars are observational estimates, 
including corrections for dust extinction: Gallego et al. (1995) 
(square); Tresse & Maddox (1998) (triangle); Gronwall (1998) 
(circle) and Sullivan et al. (2000) (cross). We have converted the 
data for the effects of differing luminosity distances and volume 
elements to correspond to the flo = 0.3, Aq = 0.7 cosmology 
assumed in our model. The solid line is the equivalent luminos- 
ity density measured in our model assuming that no photons es- 
cape (/esc, gas = 0). The dashed line is the same function for 
/esc, gas = 0.05 and the dotted line for /esc, gas = 0.2. 

namical evolution of the gas distribution around an OB as- 
sociation due to energy input by stellar winds and supernova 
from the OB association itself. Dove, ShuU & Ferrara (2000) 
have calculated the escape of ionizing photons through a dy- 
namically evolving superbubble, which is driven by an OB 
association at its centre. They find that the resulting es- 
cape fractions are slightly lower than those obtained from 
the DS94 model (since the superbubble shell is able to effec- 
tively trap radiation). Numerical solutions of the radiative 
transfer equations in disk galaxies give results in excellent 
agreement with the Stromgren sphere approach of DS94 for 
OB associations at the bright end of the luminosity function, 
but give somewhat lower escape fractions for the faintest OB 
associations, the two approaches differing by around 25% for 
a single OB star (Wood & Loeb 1999). 

Our second physical model for /esc, gas is based on De- 
vriendt et al. (1998), hereafter DSGN98. In this case, the 
ionizing stars are assumed to be uniformly mixed with the 
gas in the galaxy, and the gas is assumed to remain neutral. 
DSGN98 give an approximate analytic expression for the 
escape fraction in this case, but we have instead calculated 
the escape fraction exactly by numerical integration, for a 
specific choice for the gas density profile. 

We give details of the calculation of /esc, gas in the DS94 
and DSGN98 models in Appendix Both models contain 
one free parameter, h^/r disk, the ratio of disk scale height 
to radial scale length. We will consider the effects of vary- 
ing this parameter in §p. For starbursts, we calculate the 
escape fraction based on a simple spherical geometry, as is 
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Figure 6. The escape fraction, /esc, at 2; = as a function of 
halo mass. Thick hnes with soUd symbols show the escape frac- 
tion ignoring the effects of dust, whilst thin lines with open sym- 
bols include absorption by dust. Three models for the gas escape 
fraction are plotted: fixed escape fraction of 10% (squares), DS94 
(triangles), and DSGN98 (circles). In each case, the symbols indi- 
cate the median of the distribution of escape fractions, whilst the 
errorbars indicate the 10% and 90% intervals of the distribution. 



also described in Appendix |A[ The contribution to the total 
ionizing luminosity from bursts of star formation is small 
(< 8%) at aU redshifts. 

To summarize, we will show results from three models 
for /esc, gas as standard. These are: a model in which /esc, gas is 
held constant at 0.1; the DS94 model with OB associations 
in the disk midplane; and the DSGN98 model using our 
exact calculation of the escaping fraction. We consider the 
DS94 model to be the most realistic of our three models for 
/esc, gas, but also preseut results from the other models for 
comparison. 

In Fig. 1^ we show the variation of /esc with dark halo 
mass at 2; = for the three models. The thin and thick 
lines show the escape fraction respectively with and without 
attenuation by dust. When a halo contains more than one 
galaxy, we plot the mean /esc weighted by ionizing luminos- 
ity. At a given halo mass, halos with the lowest /esc tend 
to have the highest ionizing luminosities, as both the star 
formation rate and attenuation of photons are increased in 
galaxies with large gas contents. The three models all show 
a trend for decreasing escape fraction with increasing halo 
mass up to Mhaio ~ 10^'^ Mq. For the fixed gas escape 
fraction model, the variation in /esc is due entirely to the 
effects of dust, which can therefore be seen to be negligi- 
ble in halos less massive than ~ 10^°/i~^Mq. This decrease 
in /esc due to dust is enhanced in the other two models by 
the variation in /esc, gas, which also declines with increasing 
halo meiss. For halos more massive than W^^ Mq, the 
escape fractions rise somewhat for the variable /cse,gas mod- 
els. Note that the DSGN98 model predicts a much smaller 
escape fraction than the DS94 model at all masses. 



In Fig. |7| we plot the variation in the (ionizing 
luminosity-weighted) mean disk scale length and cold gas 
mass for galaxies in our model as a function of halo mass. 
Evidently, the decline in /esc with increasing halo mass be- 
low 10^^/i"^Mq seen in Fig. ^ is due mainly to the greater 
masses of gas found in galaxies in these halos. This rapid 
change in the mass of gas present is due to the effects of feed- 
back, which efficiently ejects gas from galaxies in low mass 
halos. Above 10^^ Mq, the mass of cold gas in galaxies 
levels off and then begins to decline as cooling becomes in- 
efficient in more massive halos. This results in an escape 
fraction increasing with halo mass for the most massive ha- 
los simulated. Although galaxy sizes increase with increasing 
halo mass, thereby reducing gas densities somewhat, this ef- 
fect is not strong enough to offset the increased cold gas 
mass in these galaxies. 

We find that the DS94 and DSGN98 models applied to 
our galaxies predict escape fractions (including the effects 
of dust) for halos of mass ~ 10^^ Mq at 2 = of ~ 
20% and x 0.2% respectively. The mean DS94 and DSGN98 
luminosity-weighted escape fractions for galaxies at z = 
are lower, being ~ 6% and ~ 0.1% respectively. However, we 
expect some variation in these values with redshift due to the 
evolution of the galaxy population. In fact, we find a rapid 
decline in both cold gas content and galaxy disk size with 
increasing redshift. In Fig. ^ we show the evolution of the 
(ionizing luminosity-weighted) mean /esc between redshifts 
and 45. All models show an initial rapid decline in /esc 
with increasing z. After this, in the constant /esc, gas model, 
the mean escape fraction increases with redshift since the 
dust content of galaxies was lower in the past. The DS94 
model shows a very gradually rising escape fraction, whilst 
the DSGN98 model has a more rapid decline. 

In our model, the contribution of stellar sources to the 
UV background is dominated by galaxies at low redshifts 
{z<,l). We find that immediately shortwards of 912A our 
DS94 model predicts a background due to stellar sources 
which is very close to that expected from QSOs (Haardt & 
Madau 1996), after including the effects of attenuation by 
the intervening IGM (Madau 1995). At shorter wavelengths 
the QSO contribution soon becomes dominant. Thus, at 
z = the combined background due to stars (from our DS94 
model) and QSOs (from Haardt & Madau 1996) is J912A ~ 
4 X 10~^^ ergs/s/cm^/Hz/ster. This is consistent with the 
upper limit of Jgi2A = 8 x 10"'^'' ergs/s/cm^/Hz/ster found 
by Vogel et al. (1995), who searched for Ho emission from 
intergalactic HI clouds. The contribution of galaxies to the 
local ionizing background has also been estimated by Gial- 
longo, Fontana & Madau (1997), based on the luminosity 
function of galaxies observed in the Canada-France Red- 
shift Survey. They estimated the galactic contribution as 
Jgj^jA ~ 5 X 10"'^'' ergs/s/cm^/Hz/ster at z = 0, assum- 
ing an escape fraction of /esc = 0.15. If we assume the 
same /esc in our model, we obtain J912A = 5.2 x 10"'^'' 
ergs/s/cm^/Hz/ster, in excellent agreement with their re- 
sult. 
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Figure 7. The variation in (ionizing luminosity-weighted) mean 
disk scale length (squares) and gas mass (triangles) for galaxies 
in our model as a function of halo mass. The symbols show the 
medians of the distributions, while the errorbars indicate the 10% 
and 90% intervals. 

4 THE FILLING FACTOR AND THE 
EVOLUTION OF THE IONIZATION 
FRONTS 

We define the filling factor, f fln, as the fraction of hydro- 
gen in the IGM (as defined in § |2.2[ ) which has been ionized. 
This is the natural quantity which serves as an indication 
of the amount of reionization in the IGM. We calculate the 
growth of the ionized region around each halo, using the 
ionizing luminosities predicted by the semi-analytic model, 
and then sum over all halos to find -Fan. We make two sim- 
plifying assumptions: (1) the radiation from each halo is 
emitted isotropically, and (2) the distribution of hydrogen is 
uniform on the scale of the ionization front and larger (but 
with small-scale clumping) . It follows that each halo by itself 
would produce a spherical ionization front. 

The mass of hydrogen ionized within the ionization 
front, M, in spherical symmetry is given by (Shapiro & 
Giroux 1987; Haiman & Loeb 1996) 

1 dM (2) _3. M 

—^b(t) — a\,'a /ciump^iH , (4) 

mn at ^ ' n ^ jj^^^ 

where nn is the comoving mean number density of hydrogen 
atoms (total, Hi and Hll) in the IGM, a{t) is the scale factor 
of the universe normalized to unity at z = 0, t is time and 
S{t) is the rate at which ionizing photons are being emitted. 
The factor /ciump(t), defined by 

/clump = (Pigm)/Pigmi (5) 

is the clumping factor for the ionized gas in the IGM (here 
piGM is the density of IGM gas at any point and pigm is the 
mean density of the IGM). Small-scale dumpiness causes 
the total recombination rate to be larger than for a uniform 
medium of the same mean density. 

The value of /ciump(i) for the ionized gas in the IGM is 
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Figure 8. The (ionizing luminosity-weighted) mean escape frac- 
tion for all galaxies as a function of redshift. Thick lines show the 
escape fraction ignoring the effects of dust, whilst thin lines in- 
clude dust. Three models for the gas escape fraction are plotted: 
constant gas escape fraction of 10% (solid line), DS94 (dotted 
hue) and DSGN98 (dashed line). 



complicated to calculate analytically. We remind the reader 
that in our picture, the IGM consists of all gas which has 
not been coUisionally ionized in halos nor become part of a 
galaxy. For a uniform IGM /dump = 1 by definition. If low 
density regions of the IGM are ionized before high density 
regions, as suggested by Miralda-Escude, Haehnelt & Rees 
(2000), then this would be similar to having /dump < 1 in 
eqn. (M), but for most purposes, /dump = 1 can be consid- 
ered as an approximate lower bound. We make two different 
estimates of the possible effects of clumping. 

For our first estimate, which we call ,fdTm™'^°'' 
sume that the photoionized gas basically traces the dark 
matter, except that gas pressure prevents it from falling 
into dark matter halos with virial temperatures smaller than 
10*K (the approximate temperature of the photo-ionized 
gas). Thus, we calculate the clumping factor as /d^mp"'^'^' ~ 
{1 + a^), where is the variance of the dark matter den- 
sity field in spheres of radius equal to the virial radius of a 
lO^K halo, is calculated from the non-linear dark matter 
power spectrum, estimated using the procedure of Peacock 
& Dodds (1996), and smoothed using a top-hat filter in real 
space. 

For our second estimate, which we call .fdu^p'j we in- 
clude the effects of coUisional ionization in halos and of re- 
moval of gas by cooling into galaxies in a way consistent 
with our definition of /igm given in equation (^). The dif- 
fuse gas in halos with virial temperatures above ~ lO^K is 
assumed to have the density profile of an isothermal sphere 
with a constant density core. The gas originally associated 
with smaller halos is assumed to be pushed out of these ha- 
los by gas pressure following photoionization, and to be in 
a uniform density component occupying the remaining vol- 
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ume. As shown in Appendix |B|, the clumping factor is then 

/int Avir 



/(halos) 
clump 



J m, smooth 



+ 



/v,smooth/K3M /k3M 



((1 - /gaOO 



x(l - a;H) 



dn 



-dMh 



(6) 



Pcf^O dMhalo 

where Mj is the mass of a halo which just retains reionized 
gas, /m, smooth IS the fractiou of the total baryonic mass in 
the uniform component, and /v, smooth is the fraction of the 
volume of the universe occupied by this gas. Here, /gai is 
the fraction of the baryonic mass in a halo in the form of 
galaxies, xh is the fraction of hydrogen in the diffuse halo gas 
which is coUisionally ionized (as in eqn. and {) indicates 
an average over all halos of mass A^haio ■ The factor ^int is 
a parameter depending only on the ratio of the size of the 
core in the gas density profile to the halo virial radius. For 
a core radius equal to one-tenth of the virial radius /int = 
3.14 (see Appendix^). This estimate of the clumping factor 
ignores the possibility of gas in the centres of halos (but not 
part of a galaxy) becoming self-shielded from the ionizing 
radiation. Such gas would not become photoionized, and so 
would not contribute to the recombination rate, resulting in 
/clump being lower than estimated here. A detailed treatment 
of the ionization and temperature structure of gas inside 
halos is beyond the scope of this work. 

Of course, before reionization, when the gas is typically 
much cooler than lO^K, gas will fall into dark matter halos 
with virial temperatures below lO^K. If all gas were in halos, 
then we would find 



/"clump 
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IGM Jo 

dn 



((1-,/gai) )(l-a:H) — TT- 

Pcilo 



dMh 



dMhalo 



(7) 



Once reionized, some of the gas in these small halos will 
flow back out of the halo as the gravitational potential is no 
longer deep enough to confine the gas. 

We consider /^'j'jj'^p' as our best estimate of the IGM 
clumping factor, at least for lower redshifts, when a signif- 
icant fraction of the gas is in halos with M > Mj, while 
r(vaiiancc) ^^^^ -j-j^g naturc of au upper limit. 

.'clump 

The clumping factors /ii''J5i°p'' and fj^^p""""^ are plotted 
as functions of redshift in Fig. bl They show fairly simi- 
lar behaviour above z ~ 10. Below this redshift, /^'^^'p"'^"' 

greatly exceeds /eiJi^°p', because it becomes dominated by 
gas in massive dark matter halos, which, on the other hand, 
contributes negligibly to /c'i'jj^]°p'' as it is coUisionally ion- 
ized. We also plot estimates of the clumping factor from 
two other papers: Valageas & Silk (1999) calculated the 
clumping factor of the baryons which have been unable to 
cool (the quantity they call Cn) using their own analytical 
model. They obtain values of /dump which are comparable to 
/clump ^ at z<,W, but are substantially larger at higher red- 
shifts. Gnedin & Ostriker (1997) performed hydrodynamical 
simulations in a cosmology similar to that which we consider, 
from which they measured /dump directly. They calculated 
two clumping factors: one for all baryons in their simulation, 
/clump ~ fiuimp^^y ^^'^ other for baryons in ionized re- 
gions only, /dump = /dump""', which is smaller. /i°°p|^"' is 
more relevant for our purposes, but may still overestimate 
the clumping of photo-ionized gas in the IGM, since it in- 




Figure 9. The gas clumping factor /dump a function of red- 
shift. The solid line shows f^^'^^"^^ whilst the dotted and dashed 

clump ' 

lines show the contributions to this quantity from gas inside and 
outside halos respectively. The dot-dashed line shows j-(™"a"'^'=) 

^ -J J clump 

Filled circles show the clumping factor calculated by Valageas & 
Silk (1999). Squares show the clumping factors determined from 
a simulation by Gnedin & Ostriker (1997) for all baryons (open 
squares) and baryons in ionized regions only (filled squares). 



eludes coUisionally ionized gas in galaxy halos. These clump- 
ing factors are everywhere lower than /dump = /chTmp"'^'^' ■ 
/ciumiT"' close to our estimate /du'^p' a^t the highest and 
lowest redshifts, but smaller in the intermediate range. 



4.1 Model results 

In Fig. |l^ we show the ionized filling factor of the IGM, -Fhn, 
as a function of redshift. Here we compute Ffin by summing 
the volumes of the Hll regions formed around each halo, 
weighted by the number of such halos per unit volume as 
given by the Press-Schechter theory. (Later we will use the 
halo mass function measured directly from an N-body sim- 
ulation to calculate ffln — see Fig. [14f). Fau will exceed 1 if 
more ionizing photons have been produced than are needed 
to completely reionize the universe. We show results for our 
three models for /esc, gas , and for three different assumptions 

about /dump. 

If we ignore the effects of dust, we find that the model 
with a constant escape fraction of 10% reionizes the Universe 
by z = 7.9 if /dump = 1 but only by z = 6.6 if /dump = 
/clump' - order to reionize the Universe by z = 5, escape 
fractions of 1.4% and 3.7% are needed for /dump = 1 and 
/clump = /dum°p' respectively. When we include the effects 
of dust, we find that gas escape fractions of 3.3% and 9.3% 
are needed to reionize by z = 5 for these two cases. 

If, instead of assuming a constant gas escape fraction, 
we use the more physically motivated DS94 model, we find 
reionization occurs at z = 6.1 if /dump = 1, but only at 
z = 4.5 if /dump ~ /chi^°p' (both estimates including dust). 
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Figure 10. The filling factor, -Fflu, as a function of redshift de- 
termined using three different models for the absorption of ioniz- 
ing photons by gas inside galaxies: constant gas escape fraction 
of 10% (solid lines), DS94 (dotted lines) and DSGN98 (dashed 
lines). All models include the effects of dust on the escape frac- 
tion. The three panels show the filling factors for three different 
assumptions about the IGM clumping factor, /dump^ no clumping 
(top panel), clumping due to virialised halos (middle panel) and 
clumping estimated from the Peacock & Dodds (1996) non-linear 
power spectrum (bottom panel). 



In the latter case, the ionized filling factor at z = 5 is only 
76%. If we assume OB associations are distributed as the gas 
in the DS94 model (as opposed to lying in the disk mid-plane 
as in our standard model) then a filling factor of 117% (i.e. 
full reionization) is achieved hy z — 5. The DSGN98 model, 
which predicts much lower escape fractions than the DS94 
model at all redshifts, is able to reionize only « 2% of the 
IGM by redshift 3, and full reionization never occurs even if 

/"clump — 1. 

We note that in the DS94 model, approximately 90% of 
the photons required for ionization are produced at z < 10. 
Thus our neglect of radiative feedback effects (which may 
reduce the number of ionizing photons produced at higher 
redshifts) is unlikely to seriously effect our determination of 
the reionization epoch. 

As both the DS94 and DSGN98 models predict quite 
low escape fractions, we have also considered a much more 
extreme model which simply assumes that /esc = + P), 
where (3 is the feedback efficiency as defined by eqn. (|l|). This 
toy model, which we will refer to as the "holes scenario", 
produces very high escaping fractions for galaxies with low 
circular speeds, and low escaping fractions for those with 
high circular speeds. A behaviour for /esc of this general 
form might result if photons are able to escape through holes 
in the galaxy disk which have been created by supernovae. 
Since dust would also be expected to be swept out of these 
holes we do not include any dust absorption in this model. 
The holes scenario produces very different results compared 
to our two physical models for the escape fraction. In this 
model /osc ~ 1 for z > 10, dropping to 45% by 2; = 0. 
Not surprisingly therefore, this model succeeds in satisfying 
the Gunn- Peterson constraint, reionizing the Universe by 
z = 11.7 if /clump = 1 and by 2: = 10.6 if /dump = .fii'irmp'- 
While this model is only a very crude attempt to consider 
a dynamically disturbed gas distribution in galaxy disks, 
it clearly demonstrates that such effects may be of great 
importance for studies of reionization. 

We have also computed the filling factor in our model 
using the clumping factors calculated by Gnedin & Ostriker 
(1997) and Valageas & Silk (1999) (as given in Fig. |). Of 
course, this is not strictly self-consistent, as their clumping 
factors are calculated from their own models for galaxy for- 
mation and reionization, which differ from ours. Using either 
of these with the DS94 model gives a reionization redshift 
comparable to that obtained using /dump = /iiu'^p'' ■ we find 
reionization at 2 = 3.6 using the Gnedin & Ostriker (1997) 
clumping factor, and z = 4.9 using that of Valageas & Silk 
(1999). 

In Fig. |l^ we show the total number of ionizing photons 
which have escaped into the IGM per unit comoving volume 
by redshift z, n^, divided by the total number of hydrogen 
nuclei in the IGM per unit comoving volume, tih (which is 
/iGM times the total number density of hydrogen nuclei). 
When this number reaches one, just enough photons have 
been emitted by galaxies to reionize the IGM completely 
if recombinations are unimportant. This criterion has been 
used previously to estimate when reionization may occur. 
Since our model includes the effects of recombinations in the 
IGM, we can judge how well this simpler criterion performs. 
If we ignore the effects of absorption by gas and dust on the 
number of ionizing photons escaping from galaxies, we find 
that, in this cosmology, our model achieves n^/nn = 1 by 
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Figure 11. The ratio of the total number of ionizing photons 
which have escaped into the IGM per comoving volume by red- 
shift z to the comoving number density of hydrogen nuclei in the 
IGM. The dot-dash line assumes that all ionizing photons can 
escape from their galaxies into the IGM. Three models for the 
gas escape fraction are also shown (each also includes the effects 
of dust): fixed escape fraction of 10% (solid line), DS94 (dotted 
line) and DSGN98 (dashed line). 



z ~ 12. When we account for the effects of dust and gas 
in galaxies, we find that the redshift at which n-y/nn = 1 
is significantly reduced, the exact value depending on the 
model for the escape fraction. With a fixed gas escape frac- 
tion of 10%, n-f/riH = 1 by 2: ~ 7, whilst for the DS94 model 
nry/riH = 1 is achieved at z ~ 6. In the DSGN98 model 
n-y/riH = 1 has not been achieved even hy z — 0. When 
we include recombinations in the IGM, the model with con- 
stant gas escape fraction reaches Fmi ~ 1 only hy z — 6.7 
for /cium^= 1, and by 2 = 5.1 for /dump = /iimnp', as shown 



in Fig. ^ showing how reionization is delayed. 

Note that while the cosmology considered here is sim- 
ilar to Model G of Baugh et al. (1998), the parameters of 
the semi-analytic model used are somewhat different. Specif- 
ically, Baugh et al. (1998) used a model in which feedback 
was much more effective in low mass halos than in our model, 
since they required their models to produce a B-band lumi- 
nosity function with a shallow faint end slope. As a result, 
the epoch at which n^/nn = 1 was much later in Model G 
of Baugh et al. (1998) than in our current model. 

In summary, we see that, even for a specific model of 
galaxy formation, the predicted epoch of reionization is sen- 
sitive to the uncertain values of the escape fraction /esc 
and the clumping factor /dump. If the clumping factor is 
as large as /d'Ji^p\ then in the case of a constant gas es- 
cape fraction /esc, gas, we need /esc, gas ^10% in our model 
to ionize the IGM by z = 5, if absorption by dust is in- 
cluded, and /esc, gas ^4% if dust is ignored. With the more 
physically-motivated DS94 and DSGN98 models, and the 
same clumping factor, at most 76% of the IGM is reionized 
by z = 5, which would be inconsistent with observations of 



the Gunn-Peterson effect. For the extreme case of a uniform 
IGM, reionization occurs by z = 6.1 even with the DS94 
model for /esc. Our "best estimate" is based on combining 
the DS94 model with /eiJi^°p'' for the IGM clumping factor. 
As already stated, this model narrowly fails to satisfy the 
Gunn-Peterson constraint at z = 5 (unless we assume that 
OB associations are distributed as the gas, rather than lying 
in the disk mid-plane), suggesting that additional sources of 
ionizing radiation are required at high redshift, either more 
stars than in our standard model, or non-stellar sources (e.g. 
quasars). However, given the theoretical uncertainties in /esc 
and /clump, we consider that this is not yet proven. 



5 SENSITIVITY OF RESULTS TO MODEL 
PARAMETERS 

We turn now to test the robustness of our results to varia- 
tions in the parameters of our galaxy formation model. To 
do this, we have varied key parameters of the models and 
determined the ionized hydrogen filling factor in each case. 
We consider several different models. The variant models 
which we consider are listed in Table In each case, we 
give the value of the parameter which is changed relative to 
the standard model given in Table |l]. 

Benson et al. (2000) have shown that normalising mod- 
els to the z = B-band luminosity function allows robust 
estimates of the z = galaxy correlation function to be 
made. Here we choose a similar constraint, forcing all mod- 
els to match the z = Ha luminosity function of Gallego et 
al. (1995) at Lua = 4 x 10*^h^^ ergs/s (note that at these 
luminosities the Gallego et al. luminosity function agrees, 
within the errorbars, with that of Sullivan et al. 2000). 
This is achieved by adjusting the value of the parameter 
T (which determines the fraction of brown dwarfs formed 
in the model). The z = Ha luminosity functions for all 
models considered are shown in Fig. |l^. Dotted lines show 
those models with Ha luminosity functions that are signifi- 
cantly different from that of the standard model (at either 
the bright or faint ends). 

In Table ^ we list escape fractions and filling factors in 
the variant models for the fixed and DS94 models for /esc 
for the case /dump = /ii'um°p^- Values of Fan may exceed 
unity, as, in some models, by z = 5 more ionizing photons 
have escaped into the IGM than are required to reionize the 
universe. 

The standard choice for the feedback efficiency, /3, 
makes feedback highly efficient in galaxies with low circular 



velocities. In this model /3 — /y"'"', where /v = Vdisk/Vhot. 
The fraction of cold gas which is reheated by supernovae 
after infinite time (a quantity with direct physical interpre- 
tation) is then 



/v 



R + P 



R + fv 



(8) 



Thus as Vdisk all gas is reheated and no stars are formed. 
For the modified feedback model, we adapt this form such 
that even in arbitrarily small potential wells not all the gas 
is reheated by supernovae. We choose 



1-R + P 



~R + fv" 



(9) 
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Table 2. Results of the variant models at z = 5. Models marked by a f are those which have very different Ho luminosity functions 
compared to the standard model (see Fig. The second column lists those parameters which differ from the standard model. The third 
column lists the value of T required for each model to match the Ha luminosity function. Columns 4-5 list the mean escape fraction 
(including the effects of dust) in the fixed and DS94 models. Finally, columns 6-7 list the filling factors at z = in the fixed and DS94 
models, for the case /dump = /chi'^p' • Note that filling factors will exceed 1 if more photons than required to reionize the universe have 
been produced. 



fescape(%) Fflu 

Model Parameter change(s) T Fixed (fesc gas = 0.1) DS94 Fixed DS94 



Standard 


None 


1.53 


3.9 


2.7 


1.08 


0.76 


It 


= 0.04 


2.87 


1.7 


0.8 


0.29 


0.13 


2t 


= 0.01 


0.36 


6.4 


7.6 


6.49 


7.60 


3t 


Ohot =0.5 


1.22 


2.8 


1.6 


1.54 


0.87 


4t 


cthot = 4.0 


1.66 


3.2 


2.0 


0.49 


0.32 


5t 


Vhot = 300 km/s 


0.95 


5.5 


4.6 


1.22 


1.03 


6t 


Vhot = 50 km/s 


0.91 


2.2 


1.3 


2.01 


1.13 


7t 


= 0.020 


1.74 


3.9 


3.3 


1.25 


1.07 


8 


e* = 0.005 


1.19 


3.9 


2.4 


1.02 


0.60 


9 


/df = 5.0 


1.47 


4.0 


2.9 


1.14 


0.83 


10 


/df = 0.2 


1.49 


2.6 


1.7 


1.04 


0.68 


lit 


/ellip = 0.05 


0.86 


2.8 


2.0 


1.78 


1.27 


12 


/dlip = 0.60 


1.63 


4.0 


2.9 


1.02 


0.72 


13 


/dyn = 5.0 


1.53 


4.0 


2.8 


1.08 


0.77 


14 


/dyn = 0.1 


1.53 


4.0 


2.9 


1.08 


0.77 


15t 


p = 0.04 


1.20 


3.0 


2.4 


1.01 


0.78 


16t 


p = 0.01 


1.46 


5.1 


3.3 


1.47 


0.96 


17 


R = 0.50 


1.86 


3.9 


2.7 


0.89 


0.63 


18t 


R = 0.10 


1.23 


4.1 


3.0 


1.34 


0.96 


19 


/iz = 0.5 


1.53 


4.0 


1.6 


1.08 


0.42 


20t 


Feedback: modified (a^ = 0.75) 


0.58 


2.7 


1.5 


3.25 


1.82 


21 


52 = 20.0 X 10^" photons/s 


1.53 


4.0 


3.4 


1.08 


0.91 


22 


52 = 0.5 X 10^0 photons/s 


1.53 


4.0 


1.4 


1.08 


0.36 


23 


Star bursts: Not included 


1.60 


4.0 


2.9 


1.04 


0.73 


24 


IMF: Salpeter (1955) 


1.45 


4.0 


2.8 


1.06 


0.75 


25t 


Gas profile: SIS 


1.59 


2.7 


1.7 


1.00 


0.64 


26t 


Recooling: allowed 


1.27 


3.3 


2.2 


1.50 


0.97 



where is an adjustable parameter. Now, as Vdisk ^ 
a fraction ar of gas is reheated, whilst a fraction 1 — 
forms stars. The standard feedback model is recovered when 
Or = 1. For the alternative feedback models considered here 
a value of = 0.75 is used. 

It should be noted that the variation having one of the 
greatest influences on the predicted fllling factors is that of 
Model 20, where we use the alternative form for feedback 
given above. This model produces an Ha luminosity func- 
tion with a very steep faint end slope, since feedback never 
becomes highly efficient, even in extremely small dark mat- 
ter halos. More ionizing photons are produced than in the 
standard model and higher fllling factors are achieved. 

Other models which alter the strength of feedback (i.e. 
Models 3, 4, 5 and 6) also cause large changes in the fllling 
factors. Models with weaker feedback (i.e. Models 3 and 6) 
result in larger filling factors as they allow more star forma- 
tion to occur in low mass galaxies (these models again pro- 
ducing a steep slope for the faint end of the Ho? luminosity 
function). The value of fib also has a strong influence on the 
filling factors as demonstrated by Models 1 and 2. Finally, 
in Models 21 and 22 we consider two alternative values of S2 
in the DS94 model. These values span the range of uncer- 
tainty for the maximum OB association luminosity in our 
own Galaxy (Dove & ShuU 1994). These models demon- 



strate that the filling factors predicted by the DS94 model 
are uncertain by a factor of at least 2 simply because of 
this uncertainty in the value of 5*2. There is, in fact, further 
uncertainty introduced as it is not clear if 5*2 represents a 
real cutoff in the luminosity function of OB associations, or 
merely a turn-over in that function. 

All of the models which significantly alter the predicted 
filliiig factors are amongst those marked with a f in Ta- 
ble y, indicating that such models do not reproduce well 
the z = Hq luminosity function, and can therefore be 
discarded as being unrealistic. With these models removed, 
our predictions for Fan are reasonably robust. Considering 
all the realistic models we find that for the fixed gas escape 
fraction of 10%, Fmi at z = 5 is 1.081^9 

(where the value 

indicates the filling factor in the standard model and the 
errors show the range found in the realistic variant models) . 
For the DS94 model we find Fm = 0-76to,ll (leaving out 
the models which vary 5*2 we find Fflii = 0.76+°-°l). Our 
conclusion that with the DS94 escape fractions reionization 
cannot happen by z = 5 if the clumping factor is as large as 
/clump' remains valid under all realistic parameter variations 
considered here. On the other hand, if the clumping factor 
is closer to the case of a uniform IGM, then reionization by 
z = 5 is possible in the DS94 model (but not in the DSGN98 
model) . 
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log(L„yiO'° Ir' ergs s"') 



Figure 12. The 2 = Ho luminosity functions of our variant 
models. The solid lines show the variant models, whilst the dotted 
lines indicate the variants models marked with a f in Table ^ 
Points with error bars are observational data from Gallego et al. 
(1995) (filled squares) and Sullivan et al. (2000) (open circles). 
These points include a correction for dust extinction. 

So far we have considered a single cosmology, namely 
ACDM. This choice was motivated by the work of Cole et al. 
(2000) and Benson et al. (2000), who have shown that the 
semi-analytic model is able to reproduce many features of 
the observed galaxy population for this cosmology. However, 
in order to explore the effects of cosmological parameters on 
reionization, we have also considered a tCDM cosmology, 
with f2 = 1, in which we use model parameters identical to 
those of Benson et al. (2000). We note that Benson et al. 
(2000) were unable to match the galaxy correlation function 
at jz = for this cosmology. 

We find that in the rCDM model, our basic results are 
unchanged, i.e. with physical models for the escape fraction, 
the IGM is reionized by z = 5 only if it is much less clumped 
than in our halo clumping model with /dump = f^^^^p^ ■ The 
escape fractions in this cosmology are actually somewhat 
higher than in the ACDM cosmology (due to a lower amount 
of gas and dust in galaxies). At 2 = 5 the DS94 model pre- 
dicts a mean escape fraction ~ 16%, whilst the DSGN98 
model predicts « 0.1%. However, the filling factors are sig- 
nificantly lower (for example, in the DS94 model Fan ~ 0.24 
at 2 = 5 when /^'Ji'^i^p' is used as compared to Fmi = 0.76 
in ACDM). This reflects the fact that many fewer ionizing 
photons are produced in this cosmology (due to the fact that 
a stronger feedback is required in order that the model fits 
the properties of galaxies at 2 = 0), and that less of the gas 
has become coUisionally ionized in virialised halos in rCDM 
than in ACDM. The only factor which works in favour of a 
higher filling factor in rCDM is that the clumping factor is 
somewhat lower. However, this is not enough to offset the 
two effects described above. 



6 SPATIAL DISTRIBUTION OF IONIZING 
SOURCES AND CMB FLUCTUATIONS 

6.1 Spatial distribution 

We now consider the temperature anisotropies imprinted on 
the microwave background by the IGM following reioniza- 
tion. These depend on the spatial and velocity correlations 
of the ionized gas. A fully self-consistent calculation of these 
correlations on the relevant scales would require very high 
resolution numerical simulations including both gas dynam- 
ics and radiative transfer (e.g. Abel & Haehnelt 1999). No 
such numerical simulation is yet available with the necessary 
combination of volume and resolution to calculate the sec- 
ondary CMB anisotropies on all angular scales of interest. 
Therefore in this paper, we calculate the spatial and veloc- 
ity distribution of the ionized gas in an approximate way, by 
combining our semi-analytical galaxy formation model with 
a high resolution N-body simulation of the dark matter. 

We have used the same ACDM simulation as Benson et 
al. (2000), described in detail by Jenkins et al. 1998, which 
has flo = 0.3, a cosmological constant Ao = 0.7, a Hubble 
constant of /i = 0.7 in units of lOOkm/s/Mpc, and which is 
normalised to produce the observed abundance of rich clus- 
ters at 2 ~ (Eke, Cole & Frenk 1996). Using the same 
semi-analytic model as employed here, Benson et al. (2000) 
were able to match the observed galaxy two-point correla- 
tion function at 2 = in this cosmology. The simulation has 
a box of length 141.3 Mpc and contains 256^ dark mat- 
ter particles, each of mass of 1.4 x 10^''h~^ Mq. We identify 
halos in this simulation using the friends-of-friends (FOF) 
algorithm with the standard linking length of 0.2, and then 
populate them with galaxies according to the semi-analytic 
model. We consider only groups consisting of 10 particles or 
more, and so resolve dark halos of mass 1.4 x 10^^ Mq or 
greater. Sources in halos which are unresolved in the simula- 
tions can produce a significant fraction of the total ionizing 
luminosity, according to the semi-analytic models. To cir- 
cumvent this problem, we add sources in unresolved halos 
into the simulation in one of two ways. The first method is 
to place the sources on randomly chosen dark matter parti- 
cles which do not belong to any resolved halo. An alterna- 
tive method is to place these sources completely at random 
within the simulation volume. This makes the unresolved 
sources completely unclustered and so is an interesting ex- 
treme case. As we will be forced to construct toy models to 
determine which regions of the simulation are ionized, the 
exact treatment of these unresolved halos will not be of great 
importance. The number of unresolved halos added to the 
simulation volume is determined from the Press-Schechter 
mass function, multiplied by a correction factor of 0.7 to 
make it match the low mass end of the N-body mass func- 
tion in ACDM at 2 = 3. 

In order to calculate the correlations between ion- 
ized regions that are needed to determine the tempera- 
ture anisotropies induced in the CMB, a simulation with 
at least the volume of this one is required. Unfortunately, 
with present computing resources, this excludes the possi- 
bility of an exact calculation of the shape and size of the 
ionized regions, which would require much higher resolu- 
tion, and also the inclusion of gas dynamics and radiative 
transfer. Therefore we have used five toy models to deter- 
mine which regions of the simulation are ionized, for a given 
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distribution of ionizing sources. These models cover a range 
of possibilities which is likely to bracket the true case, and 
provide an estimate of the present theoretical uncertainties. 

For each model, we divide the simulation volume into 
256^ cubic cells, resulting in a cell size of 0.55/i"^ Mpc. As 
the gas distribution is not homogeneous, the volume of gas 
ionized will depend on the density of gas in the ionized re- 
gion. We assume that the ionizing luminosity from the galax- 
ies in each halo all originates from the halo centre, and that 
the total mass M of gas ionized by each halo is the same 
as it would be for an IGM which is uniform on large scales, 
but with small-scale clumping /dump, as given by eqn. 
We add to this the mass of any coUisionally ionized gas in 
the halo. We then calculate the volume of the ionized region 
around each halo using M = numnV , where tTh is the mean 
IGM density within the volume V. We use several different 
toy models to calculate the spatial distribution of ionized 
gas in the simulations. In all cases, the total mass of hydro- 
gen ionized is assumed to be the same as for a homogeneous 
distribution with the specified clumping factor. 

Model A (Growing front model) Ionize a spher- 
ical volume around each halo with a radius equal to the 
ionization front radius for that halo assuming a large-scale 
uniform distribution of Hi. Since the Hi in the simulation 
is not uniformly distributed, and also because some spheres 
will overlap, the ionized volume will not contain the correct 
total mass of Hi. We therefore scale the radius of each sphere 
by a constant factor, /, and repeat the procedure. This pro- 
cess is repeated, with a new value of / each time, until the 
correct total mass of Hi has been ionized. 

Model B (High density model) In this model we 
ignore the positions of halos in the simulation. Instead we 
simply rank the cells in the simulation volume by their den- 
sity. We then completely ionize the gas in the densest cell. If 
this has not ionized enough Hi we ionize the second densest 
cell. This process is repeated until the correct total mass of 
Hi has been ionized. 

Model C (Low density model) As model B, but we 
begin by ionizing the least dense cell, and work our way up 
to cells of greater and greater density. This model mimics 
that of Miralda-Escude, Haehnelt & Rees (2000). 

Model D (Random spheres model) As Model A 
but the spheres are placed in the simulation entirely at ran- 
dom rather than on the dark matter halos. By comparing to 
Model A this model allows us to estimate the importance of 
the spatial clustering of dark matter halos. 

Model E (Boundary model) Ionize a spherical re- 
gion around each halo with a radius equal to the ioniza- 
tion front radius for that halo. This may ionize too much 
or not enough Hi depending on the density of gas around 
each source. We therefore begin adding or removing cells at 
random from the boundaries of the already ionized regions 
until the required mass of Hi is ionized. 

Fig. ^ shows six slices through the N-body simulation. 
The top left slice shows the density of all gas (which is as- 
sumed to trace the dark matter) , whilst the other slices show 
only the density of ionized gas. Model A shows particularly 
well the correlated nature of the ionizing sources (due to 
the fact that galaxies form in the high density regions of the 
dark matter), as the densest regions of the simulations are 
the ones which have become most highly ionized. 

In Fig. ^ we show the filling factor in the N-body sim- 




Figure 14. Filling factors as a function of redshift in the N- 
body simulation. Lines are plotted for the 'fixed' model with three 
different values of /osc,gas: 0.05 (solid line), 0.20 (dotted line) and 
1.00 (short-dashed line). 



ulation for the fixed /esc, gas model with different values of 
/esc, gas, for the case /clump — /c'ium°p' ■ The filling factors cal- 
culated from the simulation are always less (for a given value 

(see Fig. This is 



of ./esc, gas) than those calculated in § 4.1 



because the simulation contains fewer low mass dark matter 
halos than predicted by the Press- Schechter theory, hence it 
contains fewer ionizing sources. 

6.2 CMB Fluctuations 

The reionization of the IGM imprints secondary anisotropics 
on the CMB through Thomson scattering off free electrons 
(see, for example, Vishniac 1987; Knox, Scoccimarro & Do- 
delson 1998; Hu 2000). These anisotropics result from the 
spatially varying ionized fraction and from density and ve- 
locity variations in the ionized IGM. The calculation of these 
secondary effects involves correlation functions of density 
fiuctuations and velocity fields which are easily determined 
in our models. To predict the form of these fiuctuations, 
we first calculate the two-point correlations between ionized 
gas over the redshift range 3 to 18, assuming that gas in the 
IGM traces the dark matter density and velocity. To do this, 
we u se the 256^ grid of ionization fractions, Xc, described in 
.ll We determine in each grid cell the value of 



c = 



1 



«lo 



(10) 



where 5 is the dark matter overdensity in the cell, vios is the 
component of the mean dark matter velocity in the cell along 
the line of sight to a distant observer, and the averaging of 
Xc(l + 5) is over all cells in the simulation volume. The dark 
matter density and velocity in each cell are estimated by 
assigning the mass and velocity of each dark matter particle 
to the grid using a cloud-in-cell algorithm. 

We then compute the correlation function 
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Fully ionized box Model A 




Model B Model C 




Figure 13. The projected density of ionized gas in a slice through the ACDM N-body simulation at z = 3 shown as a grcyscalc image, 
with the densest regions being black. The slice shown has dimensions of 141.3 X 141.3 X 8.0h^^ Mpc. The total (i.e. ionized plus neutral) 
projected gas density is shown in the upper left hajid panel. The remaining panels show the projected density of ionized gas in Models 
A-E. 
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Figure 16. The effect of varying the escape fraction, /esc, gas, on 
the secondary CMB anisotropies. The curves shown are all com- 
puted using Model E, with unresolved halos placed on ungrouped 
particles. Gas escape fractions of 1.00, 0.70, 0.50, 0.20 and 0.05 
are shown as indicated in the figure. The redshift of reionization 
for each model is also indicated. The heavy solid line shows the 
primary CMB anisotropies. 



^«W = (C(x)C(x + r)).. (11) 

This correlation function is all that is needed to determine 
the spectrum of fluctuations imprinted in the CMB by the 
reionization process. A detailed description of how this spec- 
trum is computed is given in Appendix |c|. 

In Fig. |l^ we show the secondary CMB anisotropies 
calculated as described above. The left-hand panel shows 
the results for Model E with a fixed gas escape fraction 
of 0.10, for the cases that unresolved halos are placed ei- 
ther on ungrouped particles (solid line) or at random in 
the simulation volume (dotted line). The clumping factor 
is /clump ~ /c'ium°p'' i ^s wiU be used for all models considered 
in this section. The particular choice of /dump is not impor- 
tant for our conclusions about the CMB fluctuations, since 
we will consider different values of /esc, gas. The choice of 
placement scheme is seen to make little difference to the re- 
suhs, the two curves differing by <10% for 250 < ^ < 5000, 
with the difference growing to 40% by f = 20000. The size 
of the grid cell used in our calculation of the ionized gas 
correlation function corresponds to f « 25000 at z = 3, and 
the turnover around I — 10** is simply due to the cell size. 
As such, our method is unable to determine the form of the 
CMB fluctuations at higher i. 

The right-hand panel of Fig. |l^ shows the variations in 
our estimates of Ci which arise from using the flve Models 
A-E. Here the differences between the curves are larger, with 
Models B and C differing by a factor of « 2.5 at £ = 10*. 
The amplitude of the curves is affected by the strength of the 
correlations present in each model (e.g. the "high density" 
model is the most strongly correlated and has the highest 
amplitude, whilst the "low density" model has the weakest 



correlations and hence the lowest amplitude). However, the 
shapes of the curves are all very similar. 

Figure ^ examines the effect on the secondary 
anisotropies of varying the escape fraction /esc, gas . The trend 
is for increasing amplitude of anisotropy with increasing es- 
cape fraction (which results in a higher reionization red- 
shift). If, however, we boost the number of photons produced 
by increasing /esc, gas above 1 (this of course being unphysi- 
cal, but a simple way of examining the effects of producing 
more ionizing photons), little further increase in amplitude 
is seen. 

The form expected for the CMB anisotropies produced 
by patchy reionization has been calculated for a simple 
model by Gruzinov & Hu (1998). In this model, reioniza- 
tion of the universe is assumed to begin at some redshift, zi, 
and is completed (i.e. the fllling factor reaches unity) after a 
redshift interval &z. Sources are assumed to appear at ran- 
dom positions in space and to each ionize a spherical region 
of comoving radius K. Once such an ionized region has ap- 
peared it remains forever. In this model, the power l^Ct/'2-K 
is predicted to have the form of white noise at small i, since 
the ionized regions are uncorrelated. 

We compare the simple model of Gruzinov & Hu (1998) 
with our own results in Fig. |l^. Since in our model reioniza- 
tion has no well-deflned starting redshift, and ionized regions 
span a range of sizes, we simply choose values of R, Zi and 
5z in order to match the two models at the peak in the spec- 
trum (even though the position of this peak in our results 
is an artifact of our simulation resolution, we simply wish 
to demonstrate here the difference in small I slopes between 
our model and that of Gruzinov & Hu 1998). The chosen 
values of i? = 0.85/i~'^Mpc, zi = 11 and (5z = 5 are all plau- 
sible for the ionization history and sizes of ionized regions 
seen in our model (the mean comoving size of regions rang- 
ing from 1.4/i"*Mpc at z = 3 to 0.2/i"*Mpc at 2 = 18). 
Note that Knox, Scoccimarro & Dodelson (1998) calculate 
a somewhat different form for the anisotropy spectrum for 
this same model, in which the amplitude. A, is roughly half 
that found by Gruzinov & Hu (1998), and the peak in the 
spectrum occurs at slightly higher I. Despite these differ- 
ences both Gruzinov & Hu (1998) and Knox, Scoccimarro 
& Dodelson (1998) agree upon the general form of the spec- 
trum (sharp peak plus white noise at small I), and this is 
all we are interested in here. 

The Ce declines much more rapidly as ^ ^ in the 
Gruzinov & Hu (1998) model than in ours. Note that Model 
D, the random sphere model, also shows the same behaviour 
as our other models, indicating that it is not the correlated 
positions of the ionizing sources in our model which produce 
the excess power at small I. If we force all halos in our model 
to have equal ionized volumes surrounding them, whilst re- 
taining the same total fllling factor, we flnd that the excess 
power above the white noise spectrum at small I remains, so 
neither is the excess due to the range of ionizing front radii, 
R, present in our model. This excess power can therefore be 
seen to be due to the correlations in gas density and veloc- 
ity induced by gravity. In fact, if we repeat our calculations 
but ignore correlations in the gas density field (i.e. we set 
5 = everywhere) we find a CMB spectrum which has a 
slope for small I which is much closer to the Gruzinov & 
Hu (1998) white-noise slope, and which has an amplitude 
over five times lower than when density correlations are in- 
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Figure 15. The secondary CMB anisotropies generated measured from the simulation. The left hand panel shows the results for Model 
E, with a fixed gas escape fraction of 0.10. The solid line indicates a model in which unresolved halos are placed on ungrouped particles, 
whilst the dotted line shows a model with unresolved halos placed at random within the simulation volume. In the right-hand panel we 
show the results for /esc, gas = 1.00 and with unresolved halos placed on ungrouped particles. The lines show the results from the five 
different models as indicated in the figure. In each case, the heavy solid line shows the primary anisotropies in this cosmology. 




Figure 17. The effect of varying Cj, on the secondary CMB 
anisotropies. The curves shown are all computed using Model 
E, with unresolved halos placed on ungrouped particles and a gas 
escape fraction of 1.0. The solid line shows Ob = 0.02 whilst the 
dotted line shows Qb = 0.04 (Model 1 in Table |). The heavy 
solid line shows the primary CMB anisotropies. 



eluded. The remaining differences between our model and 
that of Gruzinov & Hu (1998) in this case are due to the 
correlated nature of ionized regions in our model. 

The amplitude of the secondary anisotropies also de- 
pends upon the assumed value of f2b (as this determines the 



optical depth for electron scattering) . The preferred value for 
our galaxy formation model of 57b = 0.02 is relatively low 
compared to estimates based on light element abundances 
and big bang nucleosynthesis. Fig. h/H shows the effect of in- 
creasing fib to 0.04. If the evolution of the ionized regions 
were the same in both models, we would expect the ampli- 
tude to increase by a factor of four (since it is proportional 
to the square of the baryon density). The evolution of ion- 
ized regions is actually quite similar in the two cases, and 
so the factor of four increase is seen. With this higher value 
of fib , the secondary anisotropies due to patchy reionization 
would be potentially detectable above £ ~ 3000. 

We note that a similar approach to computing the spec- 
trum of secondary anisotropies due to patchy reionization 
has been taken by Bruscoli et al. (2000). Using a differ- 
ent model of galaxy formation, Bruscoli et al. (2000) grow 
spherical ionization fronts around dark matter halos identi- 
fied in an N-body simulation, and from these they estimate 
the spectrum of CMB anisotropies produced. The simula- 
tions employed by Bruscoli et al. (2000) have higher reso- 
lution (but much smaller volume) than the GIF simulations 
used in our work. Bruscoli et al. (2000) therefore do not have 
the problem of locating unresolved halos in their simulation, 
but their calculation of secondary anisotropies is restricted 
to smaller angular scales (roughly 5 x W'^<,£<,2 x 10^) com- 
pared to ours. Furthermore, Bruscoli et al. (2000) make 
some approximations in calculating the anisotropies which 
we do not, ignoring variations in the total IGM density, and 
assuming that the ionized fraction is completely uncorre- 
lated with the velocity field. Bruscoli et al. (2000) consider 
only a single model for reionizing the simulation volume. As 
we have shown, our five toy models for the distribution of 
ionized regions lead to factors of 2-3 difference in the sec- 
ondary anisotropy amplitudes, indicating that the results 
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are not very sensitive to the model adopted for the distri- 
bution of ionized regions or to the treatment of unresolved 
halos in the simulation. Bruscoli et al. (2000) carried out 
their calculations in a different cosmology to ours, also with 
a different value of 57b, but once these differences arc taken 
into account, their results seem reasonably consistent with 
ours. 



7 DISCUSSION AND CONCLUSIONS 

We have outlined an approach to studying the reionization 
of the universe by the radiation from stars in high rcdshift 
galaxies. We have focusscd on the reionization of hydrogen, 
but the approach can be generalised to study helium reion- 
ization (e.g. Giroux & Shapiro 1994), and also to include 
radiation from quasars. Our main conclusions are: 

(i) Using a model of galaxy formation constrained by 
several observations of the local galaxy population, enough 
ionizing photons are produced to reionize the universe by 
z = 11.7. This assumes that all ionizing photons escape from 
the galaxies they originate in, and that the density of the 
IGM is uniform. Reionization is delayed until z « 10.9 in the 
case of a clumped IGM, in which gas falls into halos with 
virial temperatures exceeding lO^K. Galaxies can reionize 
such a clumped IGM by z = 5 providing that, on average, 
at least 4% of ionizing photons can escape from the galaxies 
where they arc produced. In the case of a uniform IGM, an 
escape fraction of only 1.4% is sufficient to reionize by « = 5. 

Using a physical model for the escape of ionizing radi- 
ation from galaxies, in which photons escape through "Hll 
chimneys" ionized in the gas layers in galaxy disks (Dove 
& ShuU 1994, hereafter DS94), we predict reionization by 
z = 6.1 for a uniform IGM or by « = 4.5 for a clumped 
IGM. Models which assume that all the gas in galaxy disks 
remains neutral arc unable to reionize oven a uniform IGM 
by z = 0. Using alternative estimates of the IGM clumping 
factor from Gnedin & Ostriker (1997) or Valageas & Silk 
(1999), we find reionization redshifts comparable with those 
found using our own clumping model, i.e. in the range z = 
4.5-5.0 with the DS94 model for the escape fraction. 

(ii) Once the ionizing escape fraction and IGM clump- 
ing factor have been specified, our estimates for the filling 
factor of ionized gas in the IGM axe reasonably robust, pro- 
viding that we consider only models which are successful in 
matching the Ha luminosity function of galaxies at z = 0. 
By far the greatest remaining influences on the ionized fill- 
ing factor come from the value of the baryon fraction fib and 
the prescription for feedback from supcrnovac. However, we 
have shown that altering these parameters also produces 
large changes in the z = Ha luminosity function. 

(iii) We combined our model for reionization with N- 
body simulations of the dark matter distribution in order to 
predict the spectrum of secondary anisotropics imprinted on 
the CMB by the process of reionization. The shape of this 
spectrum is almost independent of the assumptions about 
reionization, but the amplitude depends on the spatial dis- 
tribution of the ionized regions, the redshift at which reion- 
ization occurs and the baryon fraction. Wc find consider- 
ably more power in the anisotropy spectrum at small £ than 
predicted by models which do not account for the large- 
scale correlations in the gas density and velocity produced 



by gravity. Despite the uncertainty in the spatial distribu- 
tion of ionized regions, we are able to determine the ampli- 
tude of this spectrum to within a factor of three for a given 
Clh (the amplitude being proportional to Ob). The results 
found by Bruscoli ct al. (2000) using a similar technique 
are reasonably consistent with ours, once differences in Ob 
and other cosmological parameters are allowed for. 

Detection of these secondary anisotropies, which would 
constrain the reionization history of the Universe, would re- 
quire fractional temperature fluctuations of ~ 10^' to be 
measured on angular scales smaller than several arcminutes. 
Although the Planck and MAP space missions are unlikely 
to have sufBcicut sensitivity to observe such anisotropies, 
the Atacama Largo Millimeter Array is expected to be able 
to measure temperature fluctuations of the level predicted 
at £ ~ 10"* in a ten hour integration. 

Previous studies of reionization have either used an ap- 
proach similar to our own, i.e. employing some type of ana- 
lytical or semi-analytical model (e.g. Haiman & Loeb 1996; 
Valageas & Silk 1999; Chiu & Ostriker 2000; Ciardi et al. 
2000), or else have used direct hydrodynamical simulations 
(e.g. Gnedin & Ostriker 1997). While the latter technique 
can in principle follow the detailed processes of galaxy for- 
mation, gas dynamics and radiative transfer, in practice the 
resolutions attainable at present do not allow such simu- 
lations to resolve the small scales relevant to this problem. 
Furthermore, the implementation of star formation and feed- 
back in such models is far from straightforward. 

There are two main uncertainties in our approach, as in 
most others: the fraction /esc of ionizing photons that escape 
from galaxies, and the clumping factor /dump of gas in the 
IGM. Future progress depends on improving estimates of the 
effects of clumping using larger gas dynamical sinmlations, 
on better modelling of the escape of ionizing photons from 
galaxies, and on better understanding of star formation and 
supernova feedback in high redshift objects. 
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APPENDIX A: CALCULATION OF THE ESCAPING FRACTION 
Al Escaping fraction in the DS94 model: Stars in mid-plane 

In the model of Dove & ShuU (1994), hereafter DS94, ionizing photons escape from galactic disks through "Hll chimneys", 
which are holes in the neutral gas layer ionized by OB associations. The OB associations are assumed to lie in the disk 
mid-plane, and to have a distribution of ionizing luminosities dN/dS oc S"^ for Si < 5* < 5*2, (dA''/dS)dS being the number 
of associations with luminosities in the range 5* to S + dS (Kennicutt, Edgar & Hodge 1989). The gas is assumed to have a 
Gaussian vertical distribution with scaleheight h^. The fraction of Lyc photons escaping through chimneys on both sides of 
the disk at radius r is (Dove & ShuU 1994, eqn. 24) 







= < 



1 + 

where Sm is defined as 
Sni{r) = n^^'^nl exp (-2r/rdisk) h^a'^^ 
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if Si<Sm< S2 

if Sm < Si, 
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(A2) 



Here aj^ is the recombination coefficient for hydrogen for recombinations to all energy levels except the first, and we have 

(A3) 



assumed an exponential disk with radial scalelength rdisk, so that the hydrogen gas density is 
n(r, z) = no exp(-r/rdisk - z^/2h^), 



r and z being the usual cylindrical polar coordinates. 

Since 5m varies throughout the galactic disk we average the escape fraction over the entire disk, assuming that the local 
rate of star formation is proportional to the column density of the disk (Kennicutt 1989; Kennicutt 1997) and that is 
constant with radius. 

The fraction of all ionizing photons produced by the galaxy which can escape into the IGM is then given by, 
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where Sm is the value of Sm calculated for r = 0. In Figure Al we show this average escape fraction as a function of the ratio 
S^/S2 for S2/S1 = 1000. 

Our model of galaxy formation calculates the radial scale length of each galaxy's disk and also the mass of cold gas present 
in that disk, and we assume that fez /rdisk is constant. We can therefore determine no and the ratios 5^/52 and Sm/Si. Hence 
/esc, gas can be found using eqn. (A4). 



A2 Escaping fraction in DS94 model: Stars tracing gas 



In the DS94 model OB associations are assumed to lie in the midplane of the galaxy disk. If instead OB associations are 
spread throughout the gas layer, having the same vertical distribution as the cold gas, then the resulti ng e scape fraction will 
be higher than that in the DS94 model. We assume the same density profile as before, given by eqn. (A3). Consider an OB 
association emitting S ionizing photons per second, at position (r, z) in the disk. We make the assumption (as did DS94) that 
the radial variations in density can be ignored for calculating the escape fraction at radius r (which will be a valid assumption 
provided the size of the Hll region formed is much less than rdisk)- In order for any photons emitted into a cone of solid angle 
dO, which makes an angle 6 with the z-axis to escape the galaxy, the emission rate of photons into this cone must exceed the 
total recombination rate in the cone. This occurs for an ionizing luminosity 5'rcq(^), where. 



c 2 (2) f r, r \ 2 ( (z + ^cosei)^' 

Sreuid)— = noai,' exp -2 / I exp -5 — | dWiZ, 

47r V rdisk/ /n V hi 



(A5) 



which can be written as S^q{9) ~ iS^cJ/ cos"^ ^ {S^^q is the solution for cos6 > and S^cq is the solution for cos9 < 0), 
where 
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Figure Al. The average escape fraction for a galactic disk in the DS94 model with OB associations in the disk mid-plane (sohd hne) 
and distributed as the cold gas (dashed line). Both models assume S2/S1 = 1000. The dotted line (visible just above the solid line at 
small values of 8^/82) shows the effects of accounting for radial variations in the gas density in the model where OB associations are 
distributed as the gas. This line only differs noticeably from the dashed line at the lowest values of S^/S2- 



This defines two critical angles, cos et{S) = ±(5'?<;±/5')^/^ such that photons can escape the galaxy only if 6 < 9^(8) or 

e>9-{s). 

The total escaping fraction from this OB association is then given by, 
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Averaging this escape fraction over the assumed OB association luminosity function then gives a mean escape fraction of 



(A8) 







= < 



1 + 



f „Q,+ eO,+ oO,- tjO,- ~| 

i J "q _ "=i _i_ "=i _ ^r°q I 
4 \ S2 Si S2 Si J 



if S2 > Sj^e-^'-Z-d.^k 

if Si < 5S,e-2'-/'-d-k < S2 

if S'^e-^''^^''^-"'^ < Si 



(A9) 



This expression is then averaged over the galaxy disk, assuming a star formation rate proportional to the local gas density, 
to derive the mean escaping fraction for the entire galaxy. This must be done numerically 

Although we have ignored radial variations in the density of the gas when computing the escaping fraction from a single 
OB association, we find by numerical so lutio n that these variations make only a very small difference to the value of /osc,gas, 
and then only for small Sn-i/S2 (see Fig. Al). 



A3 Escaping fraction in the DSGN98 model 

In the model of Devriendt et al. (1998), DSGN98 the stars producing the Lyman continuum photons are assumed to be 
uniformly mixed with the gas in the galaxy, which is distributed in an exponential disk. All of the hydrogen in the galaxy is 
assumed to be in the form of Hi, allowing the optical depth for ionizing photons to be calculated. 



We have calculated the escaping fraction in this model exactly, using the density profile given by eqn. (A3). As in the 
case of the DS94 model with stars mixed uniformly with the gas, we begin by finding the escaping fraction as a function of 
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Figure A2. Mean escape fraction for a galaxy disk in tlie DSGN98 model as a function of tq. 



position (ro, 20) and line of sight (9, 0). For the DSGN98 model we therefore find the total optical depth in neutral hydrogen 
along the line of sight, which is 

/•oo 

T^ro,zo,6,(j)) = a-ai / n{r,z)dl, (AlO) 
Jo 

where aui ~ 6.3 x 10~^**cm^^ is the cross-section for hydrogen ionization at 912A. The quantity exp(— r) is then averaged 
over all ro, Zo, 9 and <j>, assuming that the star formation rate is proportional to the local gas density, to obtain the final 
escaping fraction. The average escape frac tion for the entire galaxy depends only on the quantity ro = nordiskCH for a given 
value of ^z/^disk, and is shown in Fig. A2 for ^z/^disk ~ 0.1. 



A4 Escaping fraction in starbursts 

In the case of a burst of star formation triggered by a major merger, we use the 83)1116 ^esCjgas as for quiescent star formation in 
the the case where /esc, gas is assumed fixed, but in the DS94 and DSGN98 models we estimate the escape fraction by assuming 
the burst has an approximately spherical geometry, throughout which star formation proceeds uniformly. We assume a sphere 
of uniform hydrogen number density, n, given by 

„ ^ 3^4as (All) 
4^»-L.stl-4mH' 

where Mgas is the mass of cold gas in the burst, rburst is the radius of the region in which the burst occurs, and mn is the mass 
of a hydrogen atom. The factor of 1.4 accounts for the presence of helium in the gas. We assume also that photons escape 
only from an outer shell of thickness I, within which the optical depth is less than 1. Therefore, 

nlam ~ 1, (A12) 

where cthi is the cross section for hydrogen ionization. The escape fraction is simply the fraction of the sphere's volume in this 
shell, i.e. 

r ^ ' burst'' /" A 1 Q\ 

^™~4V3r3_/ ^"""^ 
Substituting for I then gives 

/"esc, gas 



rburst nCTHl 

47i-rg^,,tl.4mH 



(A14) 

MgasfHl 

We take rbmst to be equal to O.lrbuige where rbuigc is the half-mass radius of the bulge formed by the merger. This choice 
is motivated by observational fact which shows that starburst activity is usually confined to the nuclear region, the size of 
which is much smaller than that of the galaxy as a whole (e.g. Sanders & Mirabel 1996 and references therein). Ricotti & 
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ShuU (1999) have carried out more elaborate calculations of escaping fractions from spherical galaxies. However, their results 
are applicable to gas in hydrostatic equilibrium with an NFW dark matter profile and so are not well suited to the case of 
starbursts. 

The star formation rate in the burst is assumed to decline exponentially, with an e-folding time equal to /dyn times the 
bulge dynamical time. Unless noted otherwise, we assume /dyn = 1 in all models. As the burst proceeds the mass of cold gas 
present, Mgas, declines as it is turned into stars. The escape fraction, /esc.gass, therefore increases during the burst, reaching 
unity as the amount of gas present drops to zero. However, as the star formation rate is declining exponentially during the 
burst only a small fraction of photons are produced whilst /esc, gas is high. 



APPENDIX B: CALCULATION OF CLUMPING FACTOR 

To estimate the clumping factor of the photoionized IGM, we make the simplifying assumption that gas in the universe can be 
split into three components — that which has fallen into dark matter halos and is coUisionally ionized or is part of a galaxy, 
that which has fallen into dark matter halos and is not coUisionally ionized, and that which has remained outside halos and 
is smoothly distributed. The first component makes no contribution to the clumping factor. We define the clumping factor as 

(Pigm) _ (Pigm) 

PiGM /iGmP 

where Pigm is the IGM gas density at any point in the universe (i.e. it does not include contributions from coUisionally ionized 
gas or galaxies), pigm ~ /igmP is the mean density of gas in the IGM and p is the mean density of all g as i n the universe 



r _ \PIGM/ _ \PIGM/ /pjix 

Jclump — _2 f2 ^2 ' \^^J 



(here /igm is the fraction of the total mass of gas in the universe which resides in the IGM, as defined in §2.2). 

Let /m, clumped be the fraction of mass in halos above the Jeans halo mass, Mj, as calculated from the Press-Schechter 
mass function for example. These halos occupy a fraction of the volume of the universe given by /v, clumped ~ /m,ciumped/Avir. 
Here, Avir is the mean density within the virial radius of a halo in units of the mean density of the Universe. The smooth 
component of gas is assumed to uniformly fill the region outside halos with M > Mj, and so has density 

Psmooth P/m, smooth/ /v, smooth; {-^^) 

where /m, smooth = 1 — /m, dumped IS the mass fraction of gas in this smooth component, and /v, smooth = 1 — /v, dumped is the 
fraction of the volume of the universe that it occupies. 

Consider next the non-coUisionally ionized gas in a single dark matter halo. Averaging over the volume of this one halo 
we obtain 

(Pdumpcd) = /int(l - /gal)^(l " a:H)^AvirP^, (B3) 

where /gai is the fraction of the baryons which have become part of galaxies within the halo, xh is the ionized fraction for 
the hydrogen in the halo gas assuming coUisional ionization equilibrium (which we take from the calculations of Sutherland 
& Dopita 1993), and /nt is a factor of order unity which depends on the shape of the halo gas density profile and is given by 

fn''" p^(r)r^dr 
Jo Pmt'" dr 

Here rvir is the virial radius of the halo, p(r) is the density profile of the difl^use gas in the halo, and pint is the mean density of 
this gas within the virial radius. We ignore any dependence of the density profile of the gas in the halo on the fraction which 
has cooled to form galaxies. Our results should be insensitive to this assumption, as /gai 1 in halos where xu is significantly 
less than unity. 

To find the contribution of gas in halos to the clumping factor, we integrate the above expression over all halos more 
massive than Mj, weighting by the volume for each halo. Adding the contribution from the smooth component, we then 
obtain 

J. /m. smooth /intAvir / / / -, j. \2\ / N2-^^halo d/T. ■, 71 ,r /T)rr\ 

/clump = -7 ■ 72 + —72 / ((l-/g-l) )(l-a:H) — ^T-TTT diV/halo, (B5) 

/v,smooth/iGM JlGM J M^ Pd lO OiWhalo 

where we have used that fact that the comoving volume of a dark matter halo of mass Mhaio is Mhaio/(AvirfioPc) (pc being 
the critical density of the universe at 2; = 0). Here ((1 — /gai)^) is averaged over all halos of mass Mhaio in our model of galaxy 
formation. 

We determine AIj by finding the mass of a dark matter halo which has a potential well deep enough that it can just hold 
onto reionized gas. This gives us the minimum mass halo within which gas collects. For the halo to just retain its gas, 

dP GMj 

d7 = 72-^('^v-), (B6) 

where rvir is the virial radius of the halo and P is the gas pressure. We approximate this as 
P GMj 



rvir 



-P(7-Vir), (B7) 



vir 
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and using the ideal gas law this becomes 

ksT GMj 47r 2 o A n , ^3 
AtniH rvir 3 

where we have used the relation Mj — 'iTvpcflo{l + z)^ Avirf'Jir/S. The virial radius is therefore 

^ ^ (1 + 

and the minimum halo mass in which gas is retained is 

471" 'J o 

We evaluate /int for the case of an isothermal profile with core radius Vc'. 

The simulations of galaxy clusters by Navarro, Frenk & White (1995) and Eke et al. (1998) show that the gas density profile 
is well described by this form. Substituting this in eqn. (B4), we find 



(B8) 



(B9) 



(BIO) 



(BU) 



arctan ■ 



arctan 



' vir ' vir 



For a typical value of Vvir/rc = 10, we therefore find /int ~ 3.14. 



(B12) 



APPENDIX C: THE SPECTRUM OF CMB SECONDARY ANISOTROPIES 

In this paper, we concentrate on the kinematic Sunyaev-Zel'dovich efi^ect which is induced by the peculiar motions (deviations 
from pure Hubble flow) of free electrons in ionized regions (Sunyaev & Zel'dovich 1980; Vishniac 1987). There exist other 
secondary sources of CMB anisotropics. However, on angular scales smaller than a few arc-minutes, the kinematic Sunyaev- 
Zel'dovich effect is likely to provide a dominant contribution. For example, it is known that the temperature anisotropies 
caused by non-linear growth of density perturbations, which are often referred to as the Rees-Sciama eff'ect or integrated 
Sachs- Wolfe effect, are of order 10~^ or less (Seljak 1996). These anisotropies depend on the the baryon bulk physical 
peculiar velocity, v, and the number density of free electrons, tIc. In our calculations of the anisotropies we assume that the 
V is equal to the bulk velocity of the dark matter and that in ionized regions is proportional to the dark matter density. 
The temperature anisotropy 6(7) — ^ observed in a given line of sight direction 7 is (e.g. Hu 2000) 

e(7,'7o) = -/ j^^-y^'^BT, (ci) 




where t] = J (1 + z)dt is conformal time with its values at recombination and present denoted, respectively, by rjrcc and 
r]o- In eqn. ( |Cl| ) we have assumed an optically thin universe. In an optically thick universe these temperature fluctuations 
are damped by a factor e~'^, where the optical depth is r = J drjaTrie/ {1 + z), where ax is the cross section for Thomson 
scattering. If the universe became instantaneously fully ionized after some redshift Zi, the relation between the optical depth 
T{rji,r]o) and Zi is approximately obtained as Zi = lOOQo {0-025/ flhh)^^^ t^^^ . Therefore, if the reionization takes place at 
z ^ lOOSlo (0.025/ flhh)^^^ , as is the case in our reionization model, then the damping factor can be neglected. 



The usual procedure to obtain the angular correlation function of temperature anisotropies in eqn. (CI) is by means of 
Limber's equation in Fourier space (see for example Peebles 1980). However, in this paper, we work in real space since we have 
the two point correlation functions of density and velocity fields directly measured in real space from N-body simulations. 

The temperature angular correlation C{9) can be written as 

rvo rvo 

C{e) = a^ drj dr7'7i7j < u'(x,7?)«^(x',7?')no(x,77)7io(x',r?') >, (C2) 

where 727" = cos^, and, x and x' refer, respectively, to comoving coordinates in the past light geodesies in the directions 7 
and 7' at r; and r;'. We write ric terms of density fluctuations 5 as 

ne(x,77) = nc{ri)xc{x,ri) [1 + (5(x,7;)] , (C3) 

where Ueiri) is the mean total (free and bound) electron number density at time rj, and 2:e(x), the ionization fraction, is unity 
in ionized regions and zero otherwise. The correlation lengths of velocity and density flelds are small compared to the Hubble 



radius so that we can approximate ndjx! ,r]') = ncijx! ,Tf) and similarly for the i;, in eqn. (S2) 
(1 + 5) 



= (1 + 5)) 



1 



«10B, (C4) 
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where «ios = JiV^ is the velocity component in the direction 7. Therefore C{6) can be written in terms of the velocity 
correlation function ^vv{y) =< uios(x)?;ios(x + y) > and the density-velocity correlation function ^cciv) =< C(x)C(x + y) >, 
both evaluated for fields at the same 77. 

Cie) = 4 r ^ r ^ne(7,)ne(77')(^o(l + S)f [Cccdx - xl) + Udx' - xl) + (C(x)z;,os(x') + C(x')^;,o.(x))] . (C5) 

The dominant contribution to C{0) is from the term involving S^q,^. The integration over yields to phase cancellation 
(Kaiser 1984; Ostriker & Vishniac 1986; Vishniac 1987). The last term in the integrand also has negligible contribution^ 
We have checked that the dominant term produces at least an order of magnitude larger anisotropics than the other ter ms. 

In flat space we use the triangle relation, jx' — xp = x'^ + t'^ — 2xx' cosd, we first carry out the integration of eqn. ( |C5[ ) 
in terms of rj' for fixed rj and 6. We compute at an average redshift z given by 1/(1 + z) = (1/(1 + zi) + 1/(1 + 22)) /2, 
which is an appropriate approximation if the correlation length is negligible relative to the horizon scale. It is straightforward 
to extend the calculation to an open geometry. 

From the temperature angular correlation C{9), we can obtain Ce as 

Ce = 2n [ dcosePi{cos6)C{e), (C6) 



where Pi{cos9) is the Legendre polynomial. 



* It is interesting that the contribution from the integral over is still dominant even if we approximate = ^ss^w, i.e. if we ignore 
any correlations between the density and velocity fields 
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